Montana Statewide nifH Analysis
Alexander Alleman
January 23, 2019
## [1] "2019-01-27-17-08-53"Working Directory
## [1] "/home/peterslab/Alex Alleman/Statewide Microbiome Analysis/Statewide analysis"set.seed(8765)Load packages
library(ggplot2)
library(data.table)
library(vegan)
library(ggvegan)
library(dplyr)
library(scales)
library(grid)
library(reshape2)
library(phyloseq)
library(ggpubr)
library(RColorBrewer)
library(ape)
library(grid)
library(knitr)
library(igraph)
library(Matrix)
library(ggnetwork)
library(intergraph)
library(parallel)
library(tinytex)Colors
farm_col<-(c("#8c510a", "#d8b365", "#f6e8c3", "#f5f5f5", "#c7eae5", "#5ab4ac", "#01665e"))
farm_col_dark<-brewer.pal(7, "Dark2")
farm_col_paired<-brewer.pal(10, "Paired")Load OTU, Taxa, and Meta data
Add OTU table with sample names on top and OTU names as row names
OTU_nifH<- read.delim(
"~/Alex Alleman/Statewide Microbiome Analysis/Statewide analysis/nifH_OTUallspring.csv",
row.names = 1)
head(OTU_nifH)[,1:10]## JZ017 JZ018 JZ019 JZ020 JZ021 JZ022 JZ023 JZ024 JZ025 JZ026
## OTU1 21 32923 48199 22 117 103 147 4054 5597 246
## OTU2 35281 5 15 1075 101 1335 26249 376 1705 256
## OTU3 11 14 13 49 20 48 9 8 16 9
## OTU4 80 4 14 2576 144 136 150 129 29299 77807
## OTU5 6 2 7 146 749 526 432 1115 13256 0
## OTU6 19317 11274 799 748 1532 348 32 21 28 2517These taxa were also created by Mr. DNA through a blast program and the NCBI database
tax_nifH<- read.delim(
"~/Alex Alleman/Statewide Microbiome Analysis/Statewide analysis/nifH_OTU_ids_2016_fixed1.txt",
row.names = 1)
head(tax_nifH)[,1:8]## kingdom phylum class order
## OTU1 k__bacteria p__firmicutes c__bacilli o__bacillales
## OTU2 k__bacteria p__firmicutes c__bacilli o__bacillales
## OTU3 k__bacteria p__proteobacteria c__alphaproteobacteria o__rhizobiales
## OTU4 k__bacteria p__firmicutes c__bacilli o__bacillales
## OTU5 k__bacteria p__proteobacteria c__alphaproteobacteria o__rhizobiales
## OTU6 k__bacteria p__proteobacteria c__alphaproteobacteria o__rhizobiales
## family genus
## OTU1 f__paenibacillaceae g__paenibacillus
## OTU2 f__paenibacillaceae g__paenibacillus
## OTU3 f__rhizobiaceae g__rhizobium
## OTU4 f__paenibacillaceae g__paenibacillus
## OTU5 f__bradyrhizobiaceae g__bradyrhizobium
## OTU6 f__bradyrhizobiaceae g__rhodopseudomonas
## species strain
## OTU1 s__paenibacillus graminis paenibacillus graminis
## OTU2 s__paenibacillus polymyxa paenibacillus polymyxa
## OTU3 s__rhizobium leguminosarum rhizobium leguminosarum
## OTU4 s__paenibacillus terrae paenibacillus terrae hpl_003
## OTU5 s__bradyrhizobium sp. bradyrhizobium sp. stm3074
## OTU6 s__rhodopseudomonas palustris rhodopseudomonas palustrisMeta data set has be placed together from all the spring and summer data with excel
meta<- read.delim(
"~/Alex Alleman/Statewide Microbiome Analysis/Statewide analysis/all_metadata_summer.csv", row.names = 1, na.strings = "NA", sep = ",",
colClasses = c(rep('factor', 7), rep('numeric', 3), rep('factor', 4), 'numeric',
rep('factor', 3), rep('numeric', 28) ) )
head(meta)[,1:5]## Site ARC Season Sample_dates Pea_variety
## JZ032 Kalispell NWARC Summer 2016-summer Delta
## JZ031 Kalispell NWARC Summer 2016-summer CDC Saffron
## JZ030 Kalispell NWARC Summer 2016-summer AC Earlystar
## JZ034 Kalispell NWARC Summer 2016-summer Majoret
## JZ033 Kalispell NWARC Summer 2016-summer DS Admiral
## JZ035 Kalispell NWARC Summer 2016-summer NavarroRemoved all Havre for analysis because samples were taken in 6" sections instead of 12"
meta2 <- meta[-c(19:48),]
sapply(meta2, class)## Site ARC Season Sample_dates
## "factor" "factor" "factor" "factor"
## Pea_variety Plot season_precip irrgation
## "factor" "factor" "numeric" "numeric"
## total_precip_irr sample_depth Date Tillage
## "numeric" "factor" "factor" "factor"
## prev_crop grain_yield elevation lat
## "factor" "numeric" "factor" "factor"
## lon Organic_Matter Moisture_Content Nitrate_Nitrite
## "factor" "numeric" "numeric" "numeric"
## Ammonia Av_Phosphorus Av_Potassium Sulfate_Sulfur
## "numeric" "numeric" "numeric" "numeric"
## pH Boron Arsenic Barium
## "numeric" "numeric" "numeric" "numeric"
## Cadmium Calcium Chromium Cobalt
## "numeric" "numeric" "numeric" "numeric"
## Copper Iron Lead Magnesium
## "numeric" "numeric" "numeric" "numeric"
## Manganese Molybdenum Nickel Phosphorus
## "numeric" "numeric" "numeric" "numeric"
## Potassium Sodium Sulfur Vanadium
## "numeric" "numeric" "numeric" "numeric"
## Zinc
## "numeric"Convert to matrix
OTU_nifH_m<-as.matrix(OTU_nifH)
tax_nifH_m<-as.matrix(tax_nifH)
meta_m<-as.matrix(meta2)
class(OTU_nifH_m)## [1] "matrix"class(tax_nifH_m)## [1] "matrix"class(meta_m)## [1] "matrix"Make phyloseq object
OTUnifH = otu_table(OTU_nifH_m, taxa_are_rows = TRUE)
TAXnifH = tax_table(tax_nifH_m)
physeq_nifH = phyloseq(OTUnifH, TAXnifH)Get physeq info
physeq_nifH## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 8821 taxa and 101 samples ]
## tax_table() Taxonomy Table: [ 8821 taxa by 8 taxonomic ranks ]Add meta data to phyloseq object
meta_phy <- sample_data(meta2)
sample_names(meta_phy)## [1] "JZ032" "JZ031" "JZ030" "JZ034" "JZ033" "JZ035" "JZ040" "JZ046"
## [9] "JZ042" "JZ044" "JZ038" "JZ036" "JZ041" "JZ047" "JZ043" "JZ045"
## [17] "JZ039" "JZ037" "JZ081" "JZ078" "JZ082" "JZ079" "JZ080" "JZ083"
## [25] "JZ105" "JZ107" "JZ103" "JZ102" "JZ106" "JZ104" "JZ084" "JZ085"
## [33] "JZ086" "JZ087" "JZ088" "JZ089" "JZ090" "JZ091" "JZ092" "JZ093"
## [41] "JZ094" "JZ095" "JZ096" "JZ097" "JZ098" "JZ099" "JZ100" "JZ101"
## [49] "JZ112" "JZ109" "JZ110" "JZ111" "JZ108" "JZ113"physeq_nifH<-merge_phyloseq(physeq_nifH, meta_phy)
physeq_nifH## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 8821 taxa and 54 samples ]
## sample_data() Sample Data: [ 54 samples by 45 sample variables ]
## tax_table() Taxonomy Table: [ 8821 taxa by 8 taxonomic ranks ]Trim data
Trim data to exclude OTUs that are not in any samples
Source of trim protocol http://evomics.org/wp-content/uploads/2016/01/phyloseq-Lab-01-Answers.html#taxa-total-counts-histogram
tdt_nifH = data.table(tax_table(physeq_nifH),
TotalCounts = taxa_sums(physeq_nifH),
OTU = taxa_names(physeq_nifH))
ggplot(tdt_nifH, aes(TotalCounts)) +
geom_histogram() +
ggtitle("Histogram of Total Counts nifH")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# How many OTUS have low count (Rare)?
tdt_nifH[(TotalCounts <= 0), .N]#zero count## [1] 1691tdt_nifH[(TotalCounts <= 1), .N]#single count## [1] 2337tdt_nifH[(TotalCounts <= 2), .N]#double count## [1] 2778# taxa cumulative sum
taxcumsum_nifH = tdt_nifH[, .N, by = TotalCounts]
setkey(taxcumsum_nifH, TotalCounts)
taxcumsum_nifH[, CumSum := cumsum(N)]
# Define the plot
pCumSum_nifH = ggplot(taxcumsum_nifH, aes(TotalCounts, CumSum)) +
geom_point() +
xlab("Filtering Threshold, Minimum Total Counts") +
ylab("OTUs Filtered") +
ggtitle("OTUs that would be filtered vs. the minimum count threshold (nifH)")
pCumSum_nifH
Zoom in
pCumSum_nifH + xlim(0, 100)## Warning: Removed 856 rows containing missing values (geom_point).
mdt_nifH = fast_melt(physeq_nifH)
prevdtnifH = mdt_nifH[, list(Prevalence = sum(count > 0),
TotalCounts = sum(count)),
by = TaxaID]ggplot(prevdtnifH, aes(Prevalence)) +
geom_histogram() +
ggtitle("Histogram of Taxa Prevalence nifH")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# How many OTUS have low prevelance (Rare)?
prevdtnifH[(Prevalence <= 0), .N]#zero ## [1] 1691prevdtnifH[(Prevalence <= 1), .N]#single ## [1] 3277prevdtnifH[(Prevalence <= 2), .N]#double ## [1] 4476prevcumsumnifH = prevdtnifH[, .N, by = Prevalence]
setkey(prevcumsumnifH, Prevalence)
prevcumsumnifH[, CumSum := cumsum(N)]
pPrevCumSumnifH = ggplot(prevcumsumnifH, aes(Prevalence, CumSum)) +
geom_point() +
xlab("Filtering Threshold, Prevalence") +
ylab("OTUs Filtered") +
ggtitle("OTUs that would be filtered vs. the minimum count threshold")
pPrevCumSumnifH
ggplot(prevdtnifH, aes(Prevalence, TotalCounts)) +
geom_point(size = 2, alpha = 0.5) +
scale_y_log10()
*** ###Trimming ####Remove less than doublets in data and prevlant in 5% of the sample ####First transform to realtive abundance
physeq_nifH_trim = filter_taxa(physeq_nifH, function(x) sum(x > 2) > (0.05*length(x)),
TRUE)
physeq_nifH_trim## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 2107 taxa and 54 samples ]
## sample_data() Sample Data: [ 54 samples by 45 sample variables ]
## tax_table() Taxonomy Table: [ 2107 taxa by 8 taxonomic ranks ]To simplfy ordination and save time we will trim the OTUs more
Remove OTUs that do not show appear more than 5 times in more than 10th of the samples
wh0_n = genefilter_sample(physeq_nifH_trim, filterfun_sample(function(x) x > 5),
A=0.1*nsamples(physeq_nifH_trim))
physeq_nifH_ord = prune_taxa(wh0_n, physeq_nifH_trim)
physeq_nifH_ord## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 804 taxa and 54 samples ]
## sample_data() Sample Data: [ 54 samples by 45 sample variables ]
## tax_table() Taxonomy Table: [ 804 taxa by 8 taxonomic ranks ]Transform to even sampling depth
physeq_nifH_ord = transform_sample_counts(physeq_nifH_ord, function(x) 1E6 * x/sum(x))
physeq_nifH_ord## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 804 taxa and 54 samples ]
## sample_data() Sample Data: [ 54 samples by 45 sample variables ]
## tax_table() Taxonomy Table: [ 804 taxa by 8 taxonomic ranks ]We have removed the majorty of the low abundance data with a remaining 804 taxa which make the data analysis much more managable.
Alpha Diveristy Analysis
Bar plots
physeq_nifH_ord_2_phylum <- tax_glom(physeq_nifH_ord, "phylum")
plot_bar(physeq_nifH_ord_2_phylum, x = "Site", fill = "phylum")+
geom_bar(aes(fill = phylum), color = "black", stat = "identity", position = "stack",
show.legend = TRUE)+
scale_fill_manual(name = "Phylum",
values=c("#7fc97f", "#beaed4", "#fdc086", "#ffff99", "#386cb0",
"#f0027f",'#e31a1c'),
#breaks=c("p_actinobacteria", "p_bacteroidetes", "p_cyanobacteria",
#"p_firmicutes", "p_proteobacteria", "p_verrucomicrobia"),
labels = c("Actinobacteria", "Bacteroidetes", "Cyanobacteria",
"Firmicutes", "Proteobacteria", "Verrucomicrobia"),
guide = guide_legend(reverse = TRUE)
)+
ggtitle("nifH Phylum Relative Abundance by Site")+
ylab("Relative Abundance (%)")+
theme(axis.text.x = element_text(angle = 90, hjust = 1), axis.text.y = element_blank(),
panel.background = element_blank())
physeq_nifH_ord_2_phylum <- tax_glom(physeq_nifH_ord, "phylum")
plot_bar(physeq_nifH_ord_2_phylum, x = "Site", fill = "phylum")+
geom_bar(aes(fill = phylum), stat = "identity", position = "stack", show.legend = TRUE)+
scale_fill_manual(name = "Phylum",
values=c("#7fc97f", "#beaed4", "#fdc086", "#ffff99",
"#386cb0", "#f0027f"),
#breaks=c("p_actinobacteria", "p_bacteroidetes", "p_cyanobacteria",
#"p_firmicutes", "p_proteobacteria", "p_verrucomicrobia"),
labels = c("Actinobacteria", "Bacteroidetes", "Cyanobacteria",
"Firmicutes", "Proteobacteria", "Verrucomicrobia"),
guide = guide_legend(reverse = TRUE)
)+
ggtitle("nifH Phylum Relative Abundance by Site")+
ylab("Relative Abundance (%)")+
theme(axis.text.x = element_text(angle = 90, hjust = 1), axis.text.y = element_blank(),
panel.background = element_blank())
Reduce phylums that contribute the less than 5% of total abundance into one group
physeq_nifH_ord_1 = transform_sample_counts(physeq_nifH_ord, function(x) x / sum(x) )
physeq_nifH_ord_phylum <- tax_glom(physeq_nifH_ord_1, "phylum")
data_nifH_phylum <- psmelt(physeq_nifH_ord_phylum)
data_nifH_phylum$phylum<-as.character(data_nifH_phylum$phylum)
data_nifH_phylum$phylum[data_nifH_phylum$Abundance<0.05]<-"<5% abdund"
#list new phylums
unique(data_nifH_phylum$phylum)## [1] "p__proteobacteria" "p__firmicutes" "p__actinobacteria"
## [4] "<5% abdund"ggplot(data = data_nifH_phylum, aes(x = Site, y = Abundance, fill = phylum))+
geom_bar(aes(fill = phylum), stat = "identity", position = "stack", show.legend = TRUE)+
scale_fill_manual(name = "Phylum",
values=c( "#33a02c", "#b2df8a","#a6cee3", "#1f78b4"),
breaks=c("p__proteobacteria", "p__firmicutes",
"p__actinobacteria","<5% abdund" ),
labels = c( "Proteobacteria", "Firmicutes",
"Actinobacteria", "<5% Abundance"),
guide = guide_legend(reverse = TRUE)
)+
ggtitle("nifH Phylum Relative Abundance by Site")+
ylab("Relative Abundance (%)")+
scale_x_discrete(labels = c("Conrad", "Corvallis", "Huntley Dryland",
"Huntley Irrigated", "Kalispell", "Mocassin",
"Sidney Dryland", "Sidney Irrigated", "Richland"))+
theme(axis.text.x = element_text(angle = 90, hjust = 1),
axis.text.y = element_blank(), panel.background = element_blank())
Publish plot to tiff
## png
## 2Different Colors
ggplot(data = data_nifH_phylum, aes(x = Site, y = Abundance, fill = phylum))+
geom_bar(aes(fill = phylum), stat = "identity", position = "stack", show.legend = TRUE)+
scale_fill_manual(name = "Phylum",
values=c( "#33a02c", "#b2df8a","#a6cee3", "#ff7f00"),
breaks=c("p__proteobacteria", "p__firmicutes",
"p__actinobacteria","<5% abdund" ),
labels = c( "Proteobacteria", "Firmicutes",
"Actinobacteria", "<5% Abundance"),
guide = guide_legend(reverse = TRUE)
)+
ggtitle("nifH Phylum Relative Abundance by Site")+
ylab("Relative Abundance")+
scale_x_discrete(labels = c("Conrad", "Corvallis",
"Huntley Dryland", "Huntley Irrigated", "Kalispell",
"Mocassin", "Sidney Dryland", "Sidney Irrigated", "Richland"))+
theme(axis.text.x = element_text(angle = 90, hjust = 1),
axis.text.y = element_blank(), panel.background = element_blank())
nifH genus barplot
physeq_nifH_ord_1 = transform_sample_counts(physeq_nifH_ord, function(x) x / sum(x) )
physeq_nifH_ord_genus <- tax_glom(physeq_nifH_ord_1, "genus")
data_nifH_genus <- psmelt(physeq_nifH_ord_genus)
data_nifH_genus$genus<-as.character(data_nifH_genus$genus)
data_nifH_genus$genus[data_nifH_genus$Abundance<0.1]<-"<10% abdund"
unique(data_nifH_genus$genus)## [1] "g__paenibacillus" "g__rhizobium" "g__geobacter"
## [4] "g__rhodopseudomonas" "g__bradyrhizobium" "g__gluconacetobacter"
## [7] "g__dechloromonas" "g__agrobacterium" "g__ideonella"
## [10] "g__azospirillum" "g__arthrobacter" "<10% abdund"ggplot(data = data_nifH_genus, aes(x = Site, y = Abundance, fill = genus))+
geom_bar(aes(fill = genus), stat = "identity", position = "stack", show.legend = TRUE)+
scale_fill_manual(name = "Genus",
values=c('#fdbf6f','#a6cee3','#b2df8a','#fb9a99','#e31a1c','#ff7f00', '#cab2d6','#6a3d9a','#ffff99','#80b1d3','#1f78b4', '#33a02c'),
breaks=c("g__paenibacillus", "g__rhizobium", "g__geobacter",
"g__rhodopseudomonas","g__bradyrhizobium",
"g__gluconacetobacter","g__dechloromonas",
"g__agrobacterium", "g__ideonella",
"g__azospirillum","g__arthrobacter", "<10% abdund"),
labels=c("Paenibacillus", "Rhizobium", "Geobacter", "Rhodopseudomonas", "Bradyrhizobium", "Gluconacetobacter", "Dechloromonas", "Agrobacterium", "Ideonella", "Azospirillum","Arthrobacter",
"<10% abdund"),
guide = guide_legend(reverse = TRUE)
)+
ggtitle("nifH Genus Relative Abundance by Site")+
ylab("Relative Abundance")+
scale_x_discrete(labels = c("Conrad", "Corvallis", "Huntley Dryland",
"Huntley Irrigated", "Kalispell", "Mocassin",
"Sidney Dryland", "Sidney Irrigated", "Richland"))+
theme(axis.text.x = element_text(angle = 90, hjust = 1), axis.text.y = element_blank(),
panel.background = element_blank())
Publish to a tiff image
## png
## 2Alpha diversity metrics
Use phyloseq internal packages to calculate the alpha diverisity
plot_richness(physeq_nifH_trim)
Simplify to just Observed and Chao1 and Shannon
plot_richness(physeq_nifH_trim, measures = c("Observed","Chao1", "Shannon"), color = "Site")
rich_nifH<-estimate_richness(physeq_nifH_trim, split = TRUE)
#First merge data sets with meta2
meta2$sample_names<-rownames(meta2)
rich_nifH$sample_names<-rownames(rich_nifH)
meta_nifH<-merge(meta2, rich_nifH, by = "sample_names")
rownames(meta_nifH)<-meta_nifH$sample_names
meta_nifH<-meta_nifH[,-1]
head(meta_nifH)[,49:54]## ACE se.ACE Shannon Simpson InvSimpson Fisher
## JZ030 1061.8638 16.02588 3.840892 0.9370300 15.880569 130.0829
## JZ031 1015.7345 15.87722 3.786340 0.9379526 16.116700 120.5734
## JZ032 985.3458 15.23833 3.241100 0.8739992 7.936459 116.9659
## JZ033 1077.3023 16.24837 3.968293 0.9463960 18.655339 133.1658
## JZ034 1044.4015 15.97121 3.962209 0.9532236 21.378282 129.1375
## JZ035 1123.6740 16.58229 3.653639 0.9231705 13.015828 133.4697Make boxplots with Observed Shannon chao1 and inverse simpson
#use ggpubr for plot
nifH_Observ<-ggboxplot(meta_nifH, x = "Site", y = "Observed",
rug = TRUE,
fill = "Site", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
nifH_Shannon<-ggboxplot(meta_nifH, x = "Site", y = "Shannon",
rug = TRUE,
fill = "Site", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
nifH_Chao<- ggboxplot(meta_nifH, x = "Site", y = "Chao1",
rug = TRUE,
fill = "Site", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
nifH_InvSim<- ggboxplot(meta_nifH, x = "Site", y = "InvSimpson",
rug = TRUE,
fill = "Site", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
alpha_nifH_fig<-ggarrange(nifH_Observ, nifH_Shannon, nifH_Chao, nifH_InvSim, ncol = 2, nrow = 2, common.legend = TRUE, legend = "bottom")
annotate_figure(alpha_nifH_fig, top = text_grob("Alpha Diversity of nifH by Site", size = 20)) 
Publish to .tiff
## png
## 2Observed- total observed OTUs
Chao1- estimate diversity and assumes that the number of observations
for a taxa has a Poisson distribution and corrects for variance
Shannon- # of OTUs (richness) scaled to the evenness
Simpson- scale of dominace probabilty of any two indviduals drawn at random beloging to the same species
Source: http://biology.kenyon.edu/courses/biol229/diversity.pdf
Just shannon plot by site
#use ggpubr for plot
ggboxplot(meta_nifH, x = "Site", y = "Shannon",
add = "mean", rug = TRUE,
fill = "Site",
title = "nifH Shannon Diversity by Site", palette = farm_col_paired, legend = "right")+
stat_compare_means(label.y = 4.3, p.adjust.method = "bonferroni")+
rotate_x_text()## Warning: Ignoring unknown parameters: p.adjust.method
Determine if distrubution is normal for each diversity test
Used the following protocol :
https://rpubs.com/dillmcfarlan/R_microbiotaSOP
#Create 2x2 plot environment so that we can see all 4 metrics at once.
par(mfrow = c(2, 2))
#Then plot each metric.
hist(rich_nifH$Shannon, main="Shannon diversity", xlab="", breaks=10)
hist(rich_nifH$InvSimpson, main="Inverse Simpson diversity", xlab="", breaks=15)
hist(rich_nifH$Chao1, main="Chao richness", xlab="", breaks=15)
hist(rich_nifH$ACE, main="ACE richness", xlab="", breaks=15)
Test for normalcy using the shapiro test. The null hypothesis for this test is that the data are normally distributed, if the p-value is greater than 0.05, then the null hypothesis is not rejected.
shapiro.test(rich_nifH$Shannon)##
## Shapiro-Wilk normality test
##
## data: rich_nifH$Shannon
## W = 0.97136, p-value = 0.2211shapiro.test(rich_nifH$InvSimpson)##
## Shapiro-Wilk normality test
##
## data: rich_nifH$InvSimpson
## W = 0.8798, p-value = 6.121e-05shapiro.test(rich_nifH$Chao1)##
## Shapiro-Wilk normality test
##
## data: rich_nifH$Chao1
## W = 0.94469, p-value = 0.01469shapiro.test(rich_nifH$ACE)##
## Shapiro-Wilk normality test
##
## data: rich_nifH$ACE
## W = 0.93376, p-value = 0.005166Use ANOVA on alpha diveristy metrics for Shannon variables becuase the other are not normal
Shannon first
aov_shannon_site_nifH <- aov(Shannon ~ Site, meta_nifH)
summary(aov_shannon_site_nifH)## Df Sum Sq Mean Sq F value Pr(>F)
## Site 8 15.624 1.9531 15.43 1.28e-10 ***
## Residuals 45 5.696 0.1266
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Site location has a signifcant effect on shannon diversity
Correct for multiple comparisons
TukeyHSD(aov_shannon_site_nifH)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Shannon ~ Site, data = meta_nifH)
##
## $Site
## diff lwr upr
## Corvallis-Conrad 0.29300868 -0.376019931 0.962037285
## Huntley Dryland -Conrad -0.50522512 -1.174253731 0.163803484
## Huntley Irrigated-Conrad 0.86418474 0.195156134 1.533213349
## Kalispell-Conrad 0.66189448 -0.007134125 1.330923091
## Mocassin-Conrad -0.44378350 -1.112812113 0.225245103
## Richland-Conrad -0.68161691 -1.350645518 -0.012588302
## Sidney Dryland-Conrad -0.66014194 -1.329170551 0.008886665
## Sidney Irrigated-Conrad 0.03893720 -0.630091409 0.707965807
## Huntley Dryland -Corvallis -0.79823380 -1.467262408 -0.129205192
## Huntley Irrigated-Corvallis 0.57117606 -0.097852543 1.240204673
## Kalispell-Corvallis 0.36888581 -0.300142802 1.037914414
## Mocassin-Corvallis -0.73679218 -1.405820790 -0.067763574
## Richland-Corvallis -0.97462559 -1.643654195 -0.305596979
## Sidney Dryland-Corvallis -0.95315062 -1.622179227 -0.284122012
## Sidney Irrigated-Corvallis -0.25407148 -0.923100086 0.414957130
## Huntley Irrigated-Huntley Dryland 1.36940986 0.700381257 2.038438473
## Kalispell-Huntley Dryland 1.16711961 0.498090998 1.836148214
## Mocassin-Huntley Dryland 0.06144162 -0.607586989 0.730470227
## Richland-Huntley Dryland -0.17639179 -0.845420395 0.492636821
## Sidney Dryland-Huntley Dryland -0.15491682 -0.823945427 0.514111789
## Sidney Irrigated-Huntley Dryland 0.54416232 -0.124866286 1.213190930
## Kalispell-Huntley Irrigated -0.20229026 -0.871318867 0.466738349
## Mocassin-Huntley Irrigated -1.30796825 -1.976996854 -0.638939638
## Richland-Huntley Irrigated -1.54580165 -2.214830260 -0.876773044
## Sidney Dryland-Huntley Irrigated -1.52432668 -2.193355292 -0.855298076
## Sidney Irrigated-Huntley Irrigated -0.82524754 -1.494276151 -0.156218935
## Mocassin-Kalispell -1.10567799 -1.774706596 -0.436649380
## Richland-Kalispell -1.34351139 -2.012540001 -0.674482785
## Sidney Dryland-Kalispell -1.32203643 -1.991065033 -0.653007818
## Sidney Irrigated-Kalispell -0.62295728 -1.291985892 0.046071324
## Richland-Mocassin -0.23783341 -0.906862013 0.431195203
## Sidney Dryland-Mocassin -0.21635844 -0.885387046 0.452670170
## Sidney Irrigated-Mocassin 0.48272070 -0.186307904 1.151749312
## Sidney Dryland-Richland 0.02147497 -0.647553640 0.690503576
## Sidney Irrigated-Richland 0.72055411 0.051525501 1.389582717
## Sidney Irrigated-Sidney Dryland 0.69907914 0.030050533 1.368107749
## p adj
## Corvallis-Conrad 0.8815206
## Huntley Dryland -Conrad 0.2786584
## Huntley Irrigated-Conrad 0.0035559
## Kalispell-Conrad 0.0545026
## Mocassin-Conrad 0.4476016
## Richland-Conrad 0.0428516
## Sidney Dryland-Conrad 0.0556618
## Sidney Irrigated-Conrad 0.9999999
## Huntley Dryland -Corvallis 0.0091627
## Huntley Irrigated-Corvallis 0.1498601
## Kalispell-Corvallis 0.6844416
## Mocassin-Corvallis 0.0211504
## Richland-Corvallis 0.0006683
## Sidney Dryland-Corvallis 0.0009315
## Sidney Irrigated-Corvallis 0.9435800
## Huntley Irrigated-Huntley Dryland 0.0000011
## Kalispell-Huntley Dryland 0.0000308
## Mocassin-Huntley Dryland 0.9999977
## Richland-Huntley Dryland 0.9940451
## Sidney Dryland-Huntley Dryland 0.9975503
## Sidney Irrigated-Huntley Dryland 0.1957237
## Kalispell-Huntley Irrigated 0.9854512
## Mocassin-Huntley Irrigated 0.0000030
## Richland-Huntley Irrigated 0.0000001
## Sidney Dryland-Huntley Irrigated 0.0000001
## Sidney Irrigated-Huntley Irrigated 0.0062523
## Mocassin-Kalispell 0.0000835
## Richland-Kalispell 0.0000017
## Sidney Dryland-Kalispell 0.0000024
## Sidney Irrigated-Kalispell 0.0858435
## Richland-Mocassin 0.9611572
## Sidney Dryland-Mocassin 0.9778887
## Sidney Irrigated-Mocassin 0.3356109
## Sidney Dryland-Richland 1.0000000
## Sidney Irrigated-Richland 0.0261606
## Sidney Irrigated-Sidney Dryland 0.0344455Plot Irrigation managment against diffrent diveristies
nifH_irr_Observ<-ggboxplot(meta_nifH, x = "Plot", y = "Observed",
rug = TRUE,
fill = "Plot", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
nifH_irr_Shannon<-ggboxplot(meta_nifH, x = "Plot", y = "Shannon",
rug = TRUE,
fill = "Plot", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
nifH_irr_Chao<- ggboxplot(meta_nifH, x = "Plot", y = "Chao1",
rug = TRUE,
fill = "Plot", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
nifH_irr_InvSim<- ggboxplot(meta_nifH, x = "Plot", y = "Simpson",
rug = TRUE,
fill = "Plot", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
alpha_nifH_irr_fig<-ggarrange(nifH_irr_Observ, nifH_irr_Shannon, nifH_irr_Chao, nifH_irr_InvSim, ncol = 2, nrow = 2, common.legend = TRUE, legend = "bottom")
annotate_figure(alpha_nifH_irr_fig, top = text_grob("Alpha Diversity of nifH by Irrigation", size = 20)) 
Looks like there is significance in the other diverisities to for irrgation (Dryland vs irrigated)
aov_observed_irr_nifH <- aov(Observed ~ Plot, meta_nifH)
summary(aov_observed_irr_nifH)## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 1351360 1351360 81.29 3.28e-12 ***
## Residuals 52 864472 16624
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Signficance in observed
aov_observed_irr_nifH <- aov(Observed ~ Plot/Site, meta_nifH)
summary(aov_observed_irr_nifH)## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 1351360 1351360 337.86 < 2e-16 ***
## Plot:Site 7 684482 97783 24.45 2.32e-13 ***
## Residuals 45 179990 4000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Plot is significant but the interaction of Plot:Site is also signficant so we cannot say there is actuall correlation
aov_Chao1_irr_nifH <- aov(Chao1 ~ Plot, meta_nifH)
summary(aov_Chao1_irr_nifH)## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 1328751 1328751 62.35 1.86e-10 ***
## Residuals 52 1108257 21313
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Signifcant in chao1
aov_Chao1_irr_nifH <- aov(Chao1 ~ Plot/Site, meta_nifH)
summary(aov_Chao1_irr_nifH)## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 1328751 1328751 234.21 < 2e-16 ***
## Plot:Site 7 852953 121850 21.48 2.11e-12 ***
## Residuals 45 255304 5673
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1aov_Simpson_irr_nifH <- aov(Simpson ~ Plot, meta_nifH)
summary(aov_Simpson_irr_nifH)## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 0.1043 0.10431 19.64 4.85e-05 ***
## Residuals 52 0.2762 0.00531
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1aov_Simpson_irr_nifH <- aov(Simpson ~ Plot/Site, meta_nifH)
summary(aov_Simpson_irr_nifH)## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 0.10431 0.10431 23.374 1.59e-05 ***
## Plot:Site 7 0.07539 0.01077 2.413 0.0347 *
## Residuals 45 0.20083 0.00446
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1All alpha diversities are signficant with irrgation method but they are nested in Site so we must look within site (Huntley and Sidney) later. It is hard to prove anything in the statewide survey with univariate ANOVAs
nifH_irr_Observ_scatter<-ggscatter(meta_nifH, x = "total_precip_irr", y = "Observed",
rug = TRUE,
add = "reg.line",
color = "Plot", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
nifH_irr_Shannon_scatter<-ggscatter(meta_nifH, x = "total_precip_irr", y = "Shannon",
rug = TRUE,
add = "reg.line",
color = "Plot", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
nifH_irr_Chao_scatter<- ggscatter(meta_nifH, x = "total_precip_irr", y = "Chao1",
rug = TRUE,
add = "reg.line",
color = "Plot", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
nifH_irr_InvSim_scatter<- ggscatter(meta_nifH, x = "total_precip_irr", y = "Simpson",
rug = TRUE,
add = "reg.line",
color = "Plot", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
alpha_nifH_irr_fig<-ggarrange(nifH_irr_Observ_scatter, nifH_irr_Shannon_scatter,
nifH_irr_Chao_scatter, nifH_irr_InvSim_scatter,
ncol = 2, nrow = 2, common.legend = TRUE, legend = "bottom")
annotate_figure(alpha_nifH_irr_fig,
top = text_grob("Alpha Diversity of nifH by Irrigation", size = 20)) 
Plot ordination
Used the following protocol from the phyloseq tutorial:
https://joey711.github.io/phyloseq/plot_ordination-examples.html
RDA Ordination
phynifH_ord_RDA <- ordinate(physeq_nifH_ord, "RDA", "bray")
plot_ordination(physeq_nifH_ord, phynifH_ord_RDA, color = "Site", shape = "Plot")+
geom_point(size = 3)+
scale_color_manual(values = farm_col_paired)+
stat_ellipse(type = "norm", linetype = 2, aes(color ="Plot"), show.legend = TRUE) +
ggtitle("nifH RDA Ordination by Site and Irrigation")+
theme_bw()
Irrigation is influencing the ordination of the principle components (PC1 is most likely comprsied of irrigation / other farm Managment)
This is really diffrenent from the irrigation RDA plot for 16s
Publish to the tiff
## png
## 2Will color with other farm managment to see if anytning is intresting.
plot_ordination(physeq_nifH_ord, phynifH_ord_RDA, color = "prev_crop", shape = "Plot")+
geom_point(size = 3)+
scale_color_manual(values = farm_col_paired)+
ggtitle("nifH RDA Ordination by Previous Crop and Irrigation")+
theme_bw()
plot_ordination(physeq_nifH_ord, phynifH_ord_RDA, color = "Tillage", shape = "Plot")+
geom_point(size = 3)+
scale_color_manual(values = farm_col_paired)+
ggtitle("nifH RDA Ordination by Previous Crop and Irrigation")+
theme_bw()
RDA plots are much diffrenent from the 16s versions but the two RDA plots from both 16s and nifH are still overlapping
Pea variety has no correlation or ordnaiton to bacterial community bray-curtis distance
plot_ordination(physeq_nifH_ord, phynifH_ord_RDA, color = "Pea_variety", shape = "Plot")+
geom_point(size = 3)+
scale_color_manual(values = farm_col_paired)+
ggtitle("nifH RDA Ordination by Previous Crop and Irrigation")+
theme_bw()
####RDA has alot of overlap
NMDS ordination
phynifH_ord_NMDS <- ordinate(physeq_nifH_ord, "NMDS", "bray")## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1618779
## Run 1 stress 0.1666296
## Run 2 stress 0.1618753
## ... New best solution
## ... Procrustes: rmse 0.001579567 max resid 0.008369045
## ... Similar to previous best
## Run 3 stress 0.1887229
## Run 4 stress 0.161878
## ... Procrustes: rmse 0.001580814 max resid 0.008390491
## ... Similar to previous best
## Run 5 stress 0.1618754
## ... Procrustes: rmse 4.945684e-05 max resid 0.0003091343
## ... Similar to previous best
## Run 6 stress 0.1618781
## ... Procrustes: rmse 0.001584558 max resid 0.008371028
## ... Similar to previous best
## Run 7 stress 0.1657018
## Run 8 stress 0.1618754
## ... Procrustes: rmse 3.746306e-05 max resid 0.0002276462
## ... Similar to previous best
## Run 9 stress 0.1887328
## Run 10 stress 0.1618754
## ... Procrustes: rmse 9.792843e-05 max resid 0.0005519838
## ... Similar to previous best
## Run 11 stress 0.1618754
## ... Procrustes: rmse 7.118828e-05 max resid 0.0004423043
## ... Similar to previous best
## Run 12 stress 0.1618779
## ... Procrustes: rmse 0.001579864 max resid 0.008364055
## ... Similar to previous best
## Run 13 stress 0.1618779
## ... Procrustes: rmse 0.00158213 max resid 0.008386758
## ... Similar to previous best
## Run 14 stress 0.1618781
## ... Procrustes: rmse 0.001578848 max resid 0.008382921
## ... Similar to previous best
## Run 15 stress 0.1618779
## ... Procrustes: rmse 0.001582478 max resid 0.008396537
## ... Similar to previous best
## Run 16 stress 0.1618754
## ... Procrustes: rmse 4.310258e-05 max resid 0.0002645992
## ... Similar to previous best
## Run 17 stress 0.1618754
## ... Procrustes: rmse 9.174953e-05 max resid 0.0004160894
## ... Similar to previous best
## Run 18 stress 0.1618754
## ... Procrustes: rmse 4.384431e-05 max resid 0.0002702027
## ... Similar to previous best
## Run 19 stress 0.1618753
## ... Procrustes: rmse 2.118693e-05 max resid 0.00010052
## ... Similar to previous best
## Run 20 stress 0.1618754
## ... Procrustes: rmse 2.938269e-05 max resid 0.0001766266
## ... Similar to previous best
## *** Solution reachedplot_ordination(physeq_nifH_ord, phynifH_ord_NMDS)
phynifH_ord_NMDS##
## Call:
## metaMDS(comm = veganifyOTU(physeq), distance = distance)
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(sqrt(veganifyOTU(physeq)))
## Distance: bray
##
## Dimensions: 2
## Stress: 0.1618753
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))'After stress test run, we get a value of 0.07 which is considered good, anything below 0.2 is acceptable.
####Plot NMDS with Site and Irrigation
plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS, color = "Site")+
geom_point(size = 3)+
scale_color_manual(values = farm_col_paired)+
ggtitle("nifH NMDS Ordination by Site")+
theme_bw()
#export physeq OTU for network analysis and saving
nifH_trim_OTU <- as(otu_table(physeq_nifH_ord), "matrix")
ID <- rownames(nifH_trim_OTU)
nifH_trim_OTU_1 <- cbind(ID, nifH_trim_OTU)
rownames(nifH_trim_OTU_1)<-c()
write.table(nifH_trim_OTU_1, file = "nifH_trim_OTU.txt", sep = "\t", row.names = FALSE)plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS, color = "Pea_variety")+
geom_point(size = 3)+
scale_color_manual(values = farm_col_paired)+
ggtitle("nifH NMDS Ordination by Pea Variety")+
theme_bw()
Pea variety does not influence the micriobiome compostion in nifH as well
plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS, shape = "Plot", color = "Site")+
geom_point(size = 3)+
scale_color_manual(values = farm_col_paired)+
ggtitle("nifH NMDS Ordination by Irrigation method")+
stat_ellipse(type = "norm", linetype = 2, aes(color = "Plot")) +
theme_bw()
Publish to tiff
## png
## 2Total water added
plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS, color = "total_precip_irr")+
geom_point(size = 3)+
scale_color_gradient(low='#05D9F6', high='#5011D1')+
ggtitle("nifH NMDS Ordination by Total Precipitation")+
theme_bw()
plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS, color = "Moisture_Content")+
geom_point(size = 3)+
scale_color_continuous(low='#05D9F6', high='#5011D1')+
ggtitle("nifH NMDS Ordination by Mositure Content")+
theme_bw()
Again total water does not ordinate well, how can we explain a diffrence between irrgation and rainfall?
Lets try some of the other variables like till, fallow, and previous crop.
plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS,color = "Tillage")+
geom_point(size = 3)+
scale_color_manual(values = farm_col_dark, breaks=c("Conventional", "Culti-roller", "No_till"),
labels=c("Conventional", "Culti Roller", "No Till"))+
ggtitle("nifH NMDS Ordination by Tillage")+
theme_bw()
Publish to .tiff
## png
## 2plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS, color = "prev_crop")+
geom_point(size = 3)+
scale_color_manual(values = c("#1B9E77", "#D95F02", "#7570B3", "#66a61E"),
name = "Previous Crop",
breaks=c("barley", "Chem_fallow", "Spring_wheat", "winter_wheat"),
labels=c("Barley", "Chemical Fallow", "Spring Wheat", "Winter Wheat"))+
ggtitle("nifH NMDS Ordination by Previous Crop")+
theme_bw()
Publish to .tiff
## png
## 2Beta dispersions
Test the diffrences in group homogeneities. Do our farm managment factors effect the homogeneitiey of the bray curtis distance?
####If a group (Site) in the MDS space are close but have diffrenent dispersion you could have a signifcant results when it is only a diffrence in dispersion.
Anderson (2006)-https://www.ncbi.nlm.nih.gov/pubmed/16542252 https://onlinelibrary.wiley.com/doi/epdf/10.1111/j.1461-0248.2006.00926.x
Irrigation beta dispersion
disp.plot <- betadisper(distance(physeq_nifH_ord, method = "bray"), meta2$Plot)
permutest(disp.plot, pairwise=TRUE, permutations=1000)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 1000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.023867 0.0238672 6.7103 1000 0.01199 *
## Residuals 52 0.184953 0.0035568
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## dryland irrigated
## dryland 0.012
## irrigated 0.012408boxplot(disp.plot)
plot(disp.plot, hull = FALSE, ellipse = TRUE)
disp.Till <- betadisper(distance(physeq_nifH_ord, method = "bray"), meta2$Tillage)
permutest(disp.Till, pairwise=TRUE, permutations=1000)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 1000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.039609 0.0198044 4.2834 1000 0.01898 *
## Residuals 51 0.235799 0.0046235
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Conventional Culti-roller No_till
## Conventional 9.9900e-04 0.0769
## Culti-roller 4.4826e-05 0.1399
## No_till 8.3245e-02 1.4099e-01boxplot(disp.Till)
plot(disp.Till, hull = FALSE, ellipse = TRUE)
disp.prev_crop <- betadisper(distance(physeq_nifH_ord, method = "bray"), meta2$prev_crop)
permutest(disp.prev_crop, pairwise=TRUE, permutations=1000)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 1000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.34683 0.115610 14.42 1000 0.000999 ***
## Residuals 50 0.40086 0.008017
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## barley Chem_fallow Spring_wheat winter_wheat
## barley 8.2917e-02 9.9900e-04 0.0010
## Chem_fallow 8.1169e-02 2.9970e-03 0.0020
## Spring_wheat 3.1442e-04 6.5930e-04 0.6963
## winter_wheat 2.7285e-05 3.6151e-05 6.7668e-01boxplot(disp.prev_crop)
plot(disp.prev_crop, hull = FALSE, ellipse = TRUE)
TukeyHSD(disp.prev_crop)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = distances ~ group, data = df)
##
## $group
## diff lwr upr p adj
## Chem_fallow-barley -0.03352468 -0.1077209 0.04067151 0.6292204
## Spring_wheat-barley -0.18749276 -0.2996669 -0.07531866 0.0002814
## winter_wheat-barley -0.22719822 -0.3393723 -0.11502412 0.0000114
## Spring_wheat-Chem_fallow -0.15396808 -0.2625802 -0.04535597 0.0023889
## winter_wheat-Chem_fallow -0.19367354 -0.3022856 -0.08506143 0.0001047
## winter_wheat-Spring_wheat -0.03970546 -0.1770901 0.09767919 0.8684302disp.site <- betadisper(distance(physeq_nifH_ord, method = "bray"), meta2$Site)
permutest(disp.site, pairwise=TRUE, permutations=1000)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 1000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 8 0.35944 0.04493 3.2652 1000 0.005994 **
## Residuals 45 0.61921 0.01376
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Conrad Corvallis Huntley Dryland Huntley Irrigated
## Conrad 0.03496503 0.03296703 0.06093906
## Corvallis 0.03646739 0.68131868 0.46953047
## Huntley Dryland 0.03003786 0.67097611 0.29470529
## Huntley Irrigated 0.06780670 0.44317758 0.29284235
## Kalispell 0.91870043 0.00041400 0.00027429 0.00246753
## Mocassin 0.52173883 0.07999231 0.06329790 0.16335615
## Richland 0.57388386 0.00439742 0.00283358 0.02313450
## Sidney Dryland 0.99250039 0.00143581 0.00099248 0.00635063
## Sidney Irrigated 0.32390060 0.21919454 0.17927660 0.38365390
## Kalispell Mocassin Richland Sidney Dryland
## Conrad 0.92007992 0.54045954 0.55844156 0.98601399
## Corvallis 0.00099900 0.08691309 0.00899101 0.00399600
## Huntley Dryland 0.00199800 0.07492507 0.00599401 0.00399600
## Huntley Irrigated 0.00899101 0.17082917 0.02897103 0.00699301
## Kalispell 0.31168831 0.26973027 0.86113886
## Mocassin 0.32083698 0.81818182 0.39160839
## Richland 0.27360130 0.82416350 0.37962038
## Sidney Dryland 0.86114730 0.40276122 0.39167738
## Sidney Irrigated 0.15341808 0.67667507 0.47279724 0.20368026
## Sidney Irrigated
## Conrad 0.3167
## Corvallis 0.2078
## Huntley Dryland 0.1738
## Huntley Irrigated 0.3956
## Kalispell 0.1409
## Mocassin 0.6643
## Richland 0.4446
## Sidney Dryland 0.2138
## Sidney Irrigatedboxplot(disp.site)
plot(disp.site, hull = FALSE, ellipse = TRUE)
TukeyHSD(disp.site)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = distances ~ group, data = df)
##
## $group
## diff lwr upr
## Corvallis-Conrad 0.1970866601 -0.02350369 0.417677013
## Huntley Dryland -Conrad 0.2056910662 -0.01489929 0.426281420
## Huntley Irrigated-Conrad 0.1730650455 -0.04752531 0.393655399
## Kalispell-Conrad -0.0092615904 -0.22985194 0.211328763
## Mocassin-Conrad 0.0683743608 -0.15221599 0.288964714
## Richland-Conrad 0.0514247120 -0.16916564 0.272015065
## Sidney Dryland-Conrad 0.0008764728 -0.21971388 0.221466826
## Sidney Irrigated-Conrad 0.1080798206 -0.11251053 0.328670174
## Huntley Dryland -Corvallis 0.0086044061 -0.21198595 0.229194759
## Huntley Irrigated-Corvallis -0.0240216145 -0.24461197 0.196568739
## Kalispell-Corvallis -0.2063482505 -0.42693860 0.014242103
## Mocassin-Corvallis -0.1287122993 -0.34930265 0.091878054
## Richland-Corvallis -0.1456619481 -0.36625230 0.074928405
## Sidney Dryland-Corvallis -0.1962101873 -0.41680054 0.024380166
## Sidney Irrigated-Corvallis -0.0890068395 -0.30959719 0.131583514
## Huntley Irrigated-Huntley Dryland -0.0326260207 -0.25321637 0.187964333
## Kalispell-Huntley Dryland -0.2149526566 -0.43554301 0.005637697
## Mocassin-Huntley Dryland -0.1373167054 -0.35790706 0.083273648
## Richland-Huntley Dryland -0.1542663542 -0.37485671 0.066323999
## Sidney Dryland-Huntley Dryland -0.2048145934 -0.42540495 0.015775760
## Sidney Irrigated-Huntley Dryland -0.0976112456 -0.31820160 0.122979108
## Kalispell-Huntley Irrigated -0.1823266359 -0.40291699 0.038263717
## Mocassin-Huntley Irrigated -0.1046906847 -0.32528104 0.115899669
## Richland-Huntley Irrigated -0.1216403336 -0.34223069 0.098950020
## Sidney Dryland-Huntley Irrigated -0.1721885728 -0.39277893 0.048401781
## Sidney Irrigated-Huntley Irrigated -0.0649852250 -0.28557558 0.155605128
## Mocassin-Kalispell 0.0776359512 -0.14295440 0.298226305
## Richland-Kalispell 0.0606863024 -0.15990405 0.281276656
## Sidney Dryland-Kalispell 0.0101380632 -0.21045229 0.230728417
## Sidney Irrigated-Kalispell 0.1173414110 -0.10324894 0.337931764
## Richland-Mocassin -0.0169496488 -0.23754000 0.203640704
## Sidney Dryland-Mocassin -0.0674978880 -0.28808824 0.153092465
## Sidney Irrigated-Mocassin 0.0397054598 -0.18088489 0.260295813
## Sidney Dryland-Richland -0.0505482392 -0.27113859 0.170042114
## Sidney Irrigated-Richland 0.0566551086 -0.16393524 0.277245462
## Sidney Irrigated-Sidney Dryland 0.1072033478 -0.11338701 0.327793701
## p adj
## Corvallis-Conrad 0.1134000
## Huntley Dryland -Conrad 0.0849905
## Huntley Irrigated-Conrad 0.2342545
## Kalispell-Conrad 1.0000000
## Mocassin-Conrad 0.9829927
## Richland-Conrad 0.9974322
## Sidney Dryland-Conrad 1.0000000
## Sidney Irrigated-Conrad 0.8019371
## Huntley Dryland -Corvallis 1.0000000
## Huntley Irrigated-Corvallis 0.9999912
## Kalispell-Corvallis 0.0830921
## Mocassin-Corvallis 0.6169037
## Richland-Corvallis 0.4537416
## Sidney Dryland-Corvallis 0.1166882
## Sidney Irrigated-Corvallis 0.9217819
## Huntley Irrigated-Huntley Dryland 0.9999076
## Kalispell-Huntley Dryland 0.0613855
## Mocassin-Huntley Dryland 0.5332525
## Richland-Huntley Dryland 0.3765001
## Sidney Dryland-Huntley Dryland 0.0875794
## Sidney Irrigated-Huntley Dryland 0.8754541
## Kalispell-Huntley Irrigated 0.1797308
## Mocassin-Huntley Irrigated 0.8276676
## Richland-Huntley Irrigated 0.6843270
## Sidney Dryland-Huntley Irrigated 0.2399571
## Sidney Irrigated-Huntley Irrigated 0.9876732
## Mocassin-Kalispell 0.9633505
## Richland-Kalispell 0.9921001
## Sidney Dryland-Kalispell 1.0000000
## Sidney Irrigated-Kalispell 0.7237181
## Richland-Mocassin 0.9999994
## Sidney Dryland-Mocassin 0.9843152
## Sidney Irrigated-Mocassin 0.9995991
## Sidney Dryland-Richland 0.9977226
## Sidney Irrigated-Richland 0.9950077
## Sidney Irrigated-Sidney Dryland 0.8087566
####The results from the betadispersion show that we have a signficant diffrence in the heterogeneity of our sites due to each of the farm amnagament factors. We see that there are sites that do not have signifcant diffrences and others that do. Beta dispersion cannot except models so we cannot detect if the diffrences are due to nestedness (prev_crop/Site).
####Do to the design of the experiment it will be hard to determine if the farm managment practices are responsible for the variation in bacterial population.
PERMANOVA (adoins)
####We will test the significance of the farm managment using the permutated-mulitvariate-ANOVA function in vegan called adonis. Adonis can test models and nestedness though the above beta dispersion test showed that most of the signficance is due to the diseprsion we cannot say much. But we can show we have the location effect.
adonis(distance(physeq_nifH_ord, method = "bray") ~Plot*prev_crop*Tillage*Pea_variety, data = meta2, permutations = 1000)##
## Call:
## adonis(formula = distance(physeq_nifH_ord, method = "bray") ~ Plot * prev_crop * Tillage * Pea_variety, data = meta2, permutations = 1000)
##
## Permutation: free
## Number of permutations: 1000
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Plot 1 1.4006 1.40059 3.8638 0.07104 0.000999 ***
## prev_crop 2 2.5341 1.26705 3.4954 0.12853 0.000999 ***
## Tillage 2 3.0720 1.53601 4.2374 0.15581 0.000999 ***
## Pea_variety 6 1.3480 0.22466 0.6198 0.06837 1.000000
## Plot:Tillage 1 1.5247 1.52471 4.2062 0.07733 0.000999 ***
## Plot:Pea_variety 6 0.9169 0.15282 0.4216 0.04651 1.000000
## prev_crop:Pea_variety 9 1.7119 0.19021 0.5247 0.08683 1.000000
## Tillage:Pea_variety 10 2.1740 0.21740 0.5997 0.11026 1.000000
## Plot:Tillage:Pea_variety 4 0.6841 0.17103 0.4718 0.03470 0.999001
## Residuals 12 4.3499 0.36249 0.22063
## Total 53 19.7163 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis(distance(physeq_nifH_ord, method = "bray") ~(Plot*prev_crop*Tillage*Pea_variety)/Site, data = meta2, permutations = 1000)##
## Call:
## adonis(formula = distance(physeq_nifH_ord, method = "bray") ~ (Plot * prev_crop * Tillage * Pea_variety)/Site, data = meta2, permutations = 1000)
##
## Permutation: free
## Number of permutations: 1000
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model
## Plot 1 1.4006 1 0
## prev_crop 2 2.5341 1 0
## Tillage 2 3.0720 2 0
## Pea_variety 6 1.3480 0 0
## Plot:Tillage 1 1.5247 2 0
## Plot:Pea_variety 6 0.9169 0 0
## prev_crop:Pea_variety 9 1.7119 0 0
## Tillage:Pea_variety 10 2.1740 0 0
## Plot:Tillage:Pea_variety 4 0.6841 0 0
## Plot:prev_crop:Tillage:Pea_variety:Site 12 4.3499 0 0
## Residuals 0 0.0000 -Inf
## Total 53 19.7163
## R2 Pr(>F)
## Plot 0.07104 1
## prev_crop 0.12853 1
## Tillage 0.15581 1
## Pea_variety 0.06837 1
## Plot:Tillage 0.07733 1
## Plot:Pea_variety 0.04651 1
## prev_crop:Pea_variety 0.08683 1
## Tillage:Pea_variety 0.11026 1
## Plot:Tillage:Pea_variety 0.03470 1
## Plot:prev_crop:Tillage:Pea_variety:Site 0.22063 1
## Residuals 0.00000
## Total 1.00000Data is nested within Site (Location effect) so the signficance in the bray-curtis dissimilarity with respect to plot is not significant with the data nested due to the lack of reproduction of conditions at each plot.
####Since there is issues with doing permutated anova over multivatrite data lets try to fit the chemical and farm managment data the NMDS orientation space using the envfit function in vegan
Model Selection
ENVFIT
Envfit does not like single variable values so we remove them
meta3<-meta2[,-c(3,4,10,11,27,29,35,38,42,46)]Will remove Site catagory liek elevation, lat, long etc that do not diffrentate between site we can call these all geographical factors as they do not change between sites.
meta3<-meta3[,-c(2,11:13)]Model fitting will be biased by chemical outliers that are in some plots the best way to avoid this is to determine the outliers (See chemical_analysis.Rmd) and remove the whole variable since functions ENVFIT and BIOENV will remove it if there are any n/a values.
####Removing Sulfate_Sulfur, Boron, Molybdenum, Potassium, Vanadium, Chromium and Sodium (From chemical_analysis.rmd)
meta3<-meta3[,-c(16,18,21,29,31)]envfitnifH <- envfit(phynifH_ord_NMDS , meta3, na.rm = TRUE, permu= 999)
envfitnifH##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## season_precip 0.30096 -0.95364 0.1693 0.008 **
## irrgation 0.54023 0.84152 0.1767 0.007 **
## total_precip_irr 0.79985 -0.60021 0.1800 0.008 **
## grain_yield -0.10337 0.99464 0.0052 0.870
## Organic_Matter -0.58271 -0.81268 0.5434 0.001 ***
## Moisture_Content 0.21634 -0.97632 0.2847 0.002 **
## Nitrate_Nitrite 0.01586 0.99987 0.0848 0.092 .
## Ammonia -0.16640 -0.98606 0.0684 0.163
## Av_Phosphorus 0.77470 0.63232 0.2592 0.001 ***
## Av_Potassium 0.98960 -0.14383 0.3455 0.001 ***
## pH -0.02160 -0.99977 0.0775 0.123
## Barium -0.09002 -0.99594 0.3593 0.001 ***
## Calcium 0.29225 -0.95634 0.2220 0.003 **
## Cobalt 0.00453 -0.99999 0.0831 0.115
## Copper 0.31097 -0.95042 0.1409 0.022 *
## Iron 0.25017 -0.96820 0.1241 0.031 *
## Magnesium 0.22694 -0.97391 0.0394 0.342
## Manganese 0.18760 -0.98225 0.1245 0.045 *
## Nickel 0.63393 -0.77339 0.1124 0.033 *
## Phosphorus 0.09644 -0.99534 0.2283 0.004 **
## Sulfur 0.10985 -0.99395 0.3851 0.001 ***
## Zinc 0.44396 -0.89605 0.2825 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## NMDS1 NMDS2
## SiteConrad 0.9037 0.2518
## SiteCorvallis 0.0861 0.7313
## SiteHuntley Dryland -0.2456 0.0338
## SiteHuntley Irrigated -0.1426 -0.0653
## SiteKalispell -1.6148 -0.2105
## SiteMocassin 0.0258 -0.7600
## SiteRichland 0.2990 -0.0971
## SiteSidney Dryland -0.1793 0.4292
## SiteSidney Irrigated 0.8678 -0.3131
## Pea_varietyAC Earlystar -0.1321 0.1464
## Pea_varietyCDC Saffron -0.0296 -0.0669
## Pea_varietyCDC Saffron 0.4876 0.1231
## Pea_varietyDelta 0.1280 -0.0638
## Pea_varietyDS Admiral 0.0304 -0.0922
## Pea_varietyMajoret -0.0978 0.0629
## Pea_varietyNavarro -0.0138 -0.0287
## Plotdryland 0.1607 -0.0285
## Plotirrigated -0.2009 0.0356
## TillageConventional -0.3088 -0.0315
## TillageCulti-roller 0.0861 0.7313
## TillageNo_till 0.1681 -0.1274
## prev_cropbarley -0.5571 0.1518
## prev_cropChem_fallow 0.1944 0.1544
## prev_cropSpring_wheat 0.8678 -0.3131
## prev_cropwinter_wheat 0.0258 -0.7600
##
## Goodness of fit:
## r2 Pr(>r)
## Site 0.8843 0.001 ***
## Pea_variety 0.0323 0.991
## Plot 0.0449 0.107
## Tillage 0.1578 0.005 **
## prev_crop 0.4008 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999The envfit function allows us to see the correlation of our envriomental vectors to the bray-curtis species dissmilarity matrix in NMDS space. This is a loose corelation to real linear corelation but it can tell us how the NMDS orinetaion is being driven.
Try a quick plot with base r and vegan for the vectors
plot(phynifH_ord_NMDS, display = "sites")
plot(envfitnifH, p.max = 0.001 )
Organic matter is resposible for 54% of the variance
Site accounts for 88% of the variance
Season precip > Irrigation
Model Selection
Bioenv
Bioenv is an iterative procedure in that links enviromental variables to community strucutre by seeking the best subset of enviromental variables that explains community strucutre
Info:
BIOENV
“The function calculates a community dissimilarity matrix using vegdist. Then it selects all possible subsets of environmental variables, scales the variables, and calculates Euclidean distances for this subset using dist. The function finds the correlation between community dissimilarities and environmental distances, and for each size of subsets, saves the best result. There are 2^p-1 subsets of p variables, and an exhaustive search may take a very, very, very long time (parameter upto offers a partial relief).”
Mantel test
“Mantel statistic is simply a correlation between entries of two dissimilarity matrices (some use cross products, but these are linearly related). However, the significance cannot be directly assessed, because there are N(N-1)/2 entries for just N observations. Mantel developed asymptotic test, but here we use permutations of N rows and columns of dissimilarity matrix.”
OTU_nifH_trim<- as(otu_table(physeq_nifH_ord ), "matrix")#Transpose the data to have sample names on rows
abund_tablenifH<-t(OTU_nifH_trim)
nrow(abund_tablenifH)## [1] 54setdiff(rownames(meta3), rownames(abund_tablenifH))## character(0)Our meta data and sample data match with 0 diffrence in rownames
Will use parallel processing to speed up calculations
#First detect amount of cores avalible
detectCores()## [1] 12#get bray-curtis distance
abund_distnifH<-vegdist(abund_tablenifH, method = "bray")
#make bioenv model against whole meta table with ".", and use gower metric to measure distance in order to incroprate the factor variables. Model will be used up to 7 vriables in final model. parelle processing on 8 cores
#nifH.bioenv <- bioenv(abund_distnifH ~ ., meta3, index = "bray", method = "pearson",
# metric = "gower", upto = 10, parallel = 10)
#summary(nifH.bioenv)#nifH.bioenvdist<-bioenvdist(nifH.bioenv, which = "best")
#mantel(nifH.bioenvdist, abund_distnifH)#nifH.bioenv$whichbest
####The best model from the inital BioEnv model selection shows that the Site explains 59% of the variation in the bray-curitis distance. But the best model explains 63% of the variation when Organic Matter and nitritie and nitrate are added in.
Try BioENV without site, this will favor farm managment variables since they are nested within.
Remove site from meta table
meta4<-meta3[,-c(1)]Rerun BioEnv
#nifH.bioenv.site.na <- bioenv(abund_distnifH ~ ., meta4, index = "bray",
# method = "pearson", metric = "gower", upto = 10, parallel = 10)
#summary(nifH.bioenv.site.na)#nifH.bioenvdist.site.na<-bioenvdist(nifH.bioenv.site.na, which = "best")
#mantel(nifH.bioenvdist.site.na, abund_distnifH)#nifH.bioenv.site.na$whichbestWhen Site is removed from the model selection we can see that there is still a bias for the Farm managment variables because they are nested with the site factor. But we do see more come out to build a bigger model. With Tillage + prev_crop + Organic_Matter + Nitrate_Nitrite + Ammonia + Av_Potassium + Manganese + Zinc.
Bioenv with chem data only
#nifH.bioenv.chem <- bioenv(abund_tablenifH ~ Organic_Matter + Moisture_Content +
# Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium +
# pH + Barium + Calcium + Cobalt + Copper + Iron + Magnesium +
# Manganese + Nickel + Phosphorus + Sulfur + Zinc, meta3,
# index = "bray", method = "pearson", metric = "gower",
# upto = 10, parallel = 10)
#summary(nifH.bioenv.chem)#nifH.bioenvdist.chem<-bioenvdist(nifH.bioenv.site.na, which = "best")
#mantel(nifH.bioenvdist.site.na, abund_distnifH)#nifH.bioenv.chem$whichbestcca/ordistep model selection
CCA model selection uses a procedure to take the a constrained distance ordination with the complete model and compare it with a unconstrainted model. The model starts with no variables then adds variables that make the best model. This is more
m1_nifH <- cca(abund_tablenifH ~ ., meta3)
m0_nifH <- cca(abund_tablenifH ~ 1, meta3)
m1_nifH## Call: cca(formula = abund_tablenifH ~ Site + Pea_variety + Plot +
## season_precip + irrgation + total_precip_irr + Tillage + prev_crop
## + grain_yield + Organic_Matter + Moisture_Content +
## Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium + pH +
## Barium + Calcium + Cobalt + Copper + Iron + Magnesium + Manganese
## + Nickel + Phosphorus + Sulfur + Zinc, data = meta3)
##
## Inertia Proportion Rank
## Total 10.7659 1.0000
## Constrained 7.8889 0.7328 33
## Unconstrained 2.8771 0.2672 20
## Inertia is scaled Chi-square
## Some constraints were aliased because they were collinear (redundant)
##
## Eigenvalues for constrained axes:
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7 CCA8 CCA9 CCA10
## 0.8070 0.6757 0.6117 0.5408 0.5291 0.4814 0.4024 0.3842 0.3332 0.3174
## CCA11 CCA12 CCA13 CCA14 CCA15 CCA16 CCA17 CCA18 CCA19 CCA20
## 0.2820 0.2391 0.2262 0.2006 0.1894 0.1707 0.1638 0.1537 0.1350 0.1302
## CCA21 CCA22 CCA23 CCA24 CCA25 CCA26 CCA27 CCA28 CCA29 CCA30
## 0.1147 0.1077 0.0963 0.0872 0.0804 0.0770 0.0695 0.0613 0.0515 0.0497
## CCA31 CCA32 CCA33
## 0.0441 0.0429 0.0330
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.29098 0.28225 0.24918 0.23022 0.20041 0.17909 0.17644 0.16025
## (Showed only 8 of all 20 unconstrained eigenvalues)m0_nifH## Call: cca(formula = abund_tablenifH ~ 1, data = meta3)
##
## Inertia Rank
## Total 10.77
## Unconstrained 10.77 53
## Inertia is scaled Chi-square
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.8189 0.6921 0.6401 0.5860 0.5584 0.5037 0.4818 0.4700
## (Showed only 8 of all 53 unconstrained eigenvalues)Ordistep
model_nifH <-ordistep(m0_nifH, scope=formula(m1_nifH))##
## Start: abund_tablenifH ~ 1
##
## Df AIC F Pr(>F)
## + Site 8 868.86 3.0678 0.005 **
## + Tillage 2 873.90 3.2418 0.005 **
## + Organic_Matter 1 874.75 3.5957 0.005 **
## + prev_crop 3 875.14 2.3857 0.005 **
## + Nitrate_Nitrite 1 875.40 2.9283 0.005 **
## + irrgation 1 875.41 2.9174 0.005 **
## + season_precip 1 875.67 2.6597 0.005 **
## + Zinc 1 875.68 2.6463 0.005 **
## + Av_Potassium 1 875.71 2.6178 0.005 **
## + total_precip_irr 1 875.74 2.5881 0.005 **
## + Nickel 1 875.88 2.4512 0.005 **
## + Cobalt 1 875.92 2.4074 0.005 **
## + Barium 1 875.93 2.3929 0.005 **
## + Iron 1 875.93 2.3923 0.005 **
## + Copper 1 875.99 2.3375 0.005 **
## + Plot 1 876.00 2.3303 0.005 **
## + Magnesium 1 876.00 2.3264 0.005 **
## + Av_Phosphorus 1 876.01 2.3202 0.005 **
## + Phosphorus 1 876.07 2.2549 0.005 **
## + Moisture_Content 1 876.18 2.1464 0.005 **
## + Manganese 1 876.18 2.1462 0.005 **
## + grain_yield 1 876.27 2.0533 0.005 **
## + Ammonia 1 876.47 1.8528 0.005 **
## + Sulfur 1 876.47 1.8514 0.005 **
## + Calcium 1 876.57 1.7584 0.005 **
## + pH 1 876.84 1.4897 0.015 *
## + Pea_variety 6 882.82 0.8462 0.980
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Site
##
## Df AIC F Pr(>F)
## - Site 8 876.36 3.0678 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Nitrate_Nitrite 1 868.08 2.3207 0.005 **
## + Ammonia 1 869.09 1.4610 0.030 *
## + Av_Phosphorus 1 869.18 1.3900 0.060 .
## + Pea_variety 6 872.61 1.0730 0.185
## + Cobalt 1 869.64 1.0051 0.370
## + Manganese 1 869.65 0.9978 0.420
## + grain_yield 1 869.68 0.9722 0.445
## + Iron 1 869.74 0.9187 0.460
## + Phosphorus 1 869.79 0.8751 0.565
## + Copper 1 869.81 0.8624 0.565
## + Zinc 1 869.77 0.8984 0.575
## + Barium 1 869.78 0.8907 0.600
## + Av_Potassium 1 869.86 0.8240 0.705
## + Magnesium 1 869.87 0.8137 0.705
## + Nickel 1 869.82 0.8508 0.730
## + pH 1 869.81 0.8615 0.750
## + Organic_Matter 1 870.05 0.6648 0.900
## + Moisture_Content 1 869.98 0.7224 0.910
## + Calcium 1 870.08 0.6382 0.925
## + Sulfur 1 870.03 0.6823 0.950
## + Plot 0 868.86
## + season_precip 0 868.86
## + irrgation 0 868.86
## + total_precip_irr 0 868.86
## + Tillage 0 868.86
## + prev_crop 0 868.86
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Site + Nitrate_Nitrite
##
## Df AIC F Pr(>F)
## - Nitrate_Nitrite 1 868.86 2.3207 0.005 **
## - Site 8 875.40 2.9709 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Ammonia 1 868.19 1.5305 0.025 *
## + Av_Phosphorus 1 868.30 1.4419 0.070 .
## + Pea_variety 6 871.61 1.0759 0.190
## + Cobalt 1 868.83 1.0065 0.445
## + Zinc 1 868.92 0.9313 0.450
## + Manganese 1 868.81 1.0230 0.455
## + Iron 1 868.90 0.9530 0.495
## + Phosphorus 1 868.90 0.9475 0.525
## + Copper 1 868.95 0.9077 0.550
## + Barium 1 868.93 0.9309 0.580
## + Av_Potassium 1 868.97 0.8946 0.600
## + Magnesium 1 869.03 0.8495 0.665
## + pH 1 868.99 0.8775 0.690
## + Nickel 1 869.05 0.8322 0.765
## + grain_yield 1 869.10 0.7851 0.840
## + Organic_Matter 1 869.20 0.7088 0.875
## + Moisture_Content 1 869.15 0.7450 0.900
## + Sulfur 1 869.21 0.6989 0.920
## + Calcium 1 869.26 0.6602 0.925
## + Plot 0 868.08
## + season_precip 0 868.08
## + irrgation 0 868.08
## + total_precip_irr 0 868.08
## + Tillage 0 868.08
## + prev_crop 0 868.08
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Site + Nitrate_Nitrite + Ammonia
##
## Df AIC F Pr(>F)
## - Ammonia 1 868.08 1.5305 0.035 *
## - Nitrate_Nitrite 1 869.09 2.3726 0.005 **
## - Site 8 875.48 2.8982 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Av_Phosphorus 1 868.52 1.3192 0.080 .
## + Pea_variety 6 871.49 1.0780 0.155
## + Iron 1 868.97 0.9609 0.425
## + Cobalt 1 868.93 0.9935 0.470
## + Manganese 1 868.94 0.9858 0.475
## + Zinc 1 869.01 0.9343 0.505
## + Barium 1 869.00 0.9394 0.515
## + Copper 1 869.03 0.9155 0.535
## + Phosphorus 1 868.98 0.9541 0.545
## + pH 1 869.03 0.9173 0.620
## + Magnesium 1 869.10 0.8613 0.670
## + Av_Potassium 1 869.11 0.8515 0.720
## + Nickel 1 869.12 0.8419 0.740
## + grain_yield 1 869.16 0.8129 0.825
## + Organic_Matter 1 869.28 0.7172 0.885
## + Sulfur 1 869.27 0.7237 0.905
## + Calcium 1 869.33 0.6769 0.935
## + Moisture_Content 1 869.30 0.6981 0.960
## + Plot 0 868.19
## + season_precip 0 868.19
## + irrgation 0 868.19
## + total_precip_irr 0 868.19
## + Tillage 0 868.19
## + prev_crop 0 868.19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cca without site
m1_nifH_site_na <- cca(abund_tablenifH ~ ., meta4)
m0_nifH_site_na <- cca(abund_tablenifH ~ 1, meta4)
m1_nifH_site_na## Call: cca(formula = abund_tablenifH ~ Pea_variety + Plot +
## season_precip + irrgation + total_precip_irr + Tillage + prev_crop
## + grain_yield + Organic_Matter + Moisture_Content +
## Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium + pH +
## Barium + Calcium + Cobalt + Copper + Iron + Magnesium + Manganese
## + Nickel + Phosphorus + Sulfur + Zinc, data = meta4)
##
## Inertia Proportion Rank
## Total 10.7659 1.0000
## Constrained 7.6960 0.7148 32
## Unconstrained 3.0700 0.2852 21
## Inertia is scaled Chi-square
## Some constraints were aliased because they were collinear (redundant)
##
## Eigenvalues for constrained axes:
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7 CCA8 CCA9 CCA10
## 0.8067 0.6725 0.6096 0.5394 0.5055 0.4714 0.3948 0.3835 0.3273 0.3114
## CCA11 CCA12 CCA13 CCA14 CCA15 CCA16 CCA17 CCA18 CCA19 CCA20
## 0.2705 0.2277 0.2256 0.1944 0.1815 0.1687 0.1585 0.1494 0.1316 0.1247
## CCA21 CCA22 CCA23 CCA24 CCA25 CCA26 CCA27 CCA28 CCA29 CCA30
## 0.1138 0.1030 0.0960 0.0871 0.0803 0.0698 0.0633 0.0572 0.0498 0.0445
## CCA31 CCA32
## 0.0430 0.0334
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.29949 0.28234 0.26627 0.23089 0.20393 0.19999 0.17647 0.17203
## (Showed only 8 of all 21 unconstrained eigenvalues)m0_nifH_site_na## Call: cca(formula = abund_tablenifH ~ 1, data = meta4)
##
## Inertia Rank
## Total 10.77
## Unconstrained 10.77 53
## Inertia is scaled Chi-square
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.8189 0.6921 0.6401 0.5860 0.5584 0.5037 0.4818 0.4700
## (Showed only 8 of all 53 unconstrained eigenvalues)Ordistep
model_nifH <-ordistep(m0_nifH_site_na, scope=formula(m1_nifH_site_na))##
## Start: abund_tablenifH ~ 1
##
## Df AIC F Pr(>F)
## + Tillage 2 873.90 3.2418 0.005 **
## + Organic_Matter 1 874.75 3.5957 0.005 **
## + prev_crop 3 875.14 2.3857 0.005 **
## + Nitrate_Nitrite 1 875.40 2.9283 0.005 **
## + irrgation 1 875.41 2.9174 0.005 **
## + season_precip 1 875.67 2.6597 0.005 **
## + Zinc 1 875.68 2.6463 0.005 **
## + Av_Potassium 1 875.71 2.6178 0.005 **
## + total_precip_irr 1 875.74 2.5881 0.005 **
## + Nickel 1 875.88 2.4512 0.005 **
## + Cobalt 1 875.92 2.4074 0.005 **
## + Barium 1 875.93 2.3929 0.005 **
## + Iron 1 875.93 2.3923 0.005 **
## + Copper 1 875.99 2.3375 0.005 **
## + Plot 1 876.00 2.3303 0.005 **
## + Magnesium 1 876.00 2.3264 0.005 **
## + Av_Phosphorus 1 876.01 2.3202 0.005 **
## + Phosphorus 1 876.07 2.2549 0.005 **
## + Moisture_Content 1 876.18 2.1464 0.005 **
## + Manganese 1 876.18 2.1462 0.005 **
## + grain_yield 1 876.27 2.0533 0.005 **
## + Ammonia 1 876.47 1.8528 0.005 **
## + Calcium 1 876.57 1.7584 0.005 **
## + Sulfur 1 876.47 1.8514 0.010 **
## + pH 1 876.84 1.4897 0.020 *
## + Pea_variety 6 882.82 0.8462 0.980
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Tillage
##
## Df AIC F Pr(>F)
## - Tillage 2 876.36 3.2418 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + prev_crop 3 871.82 2.5825 0.005 **
## + Organic_Matter 1 872.19 3.5560 0.005 **
## + total_precip_irr 1 872.60 3.1476 0.005 **
## + irrgation 1 872.87 2.8842 0.005 **
## + Av_Potassium 1 872.95 2.8096 0.005 **
## + Av_Phosphorus 1 872.98 2.7771 0.005 **
## + season_precip 1 873.16 2.6012 0.005 **
## + Magnesium 1 873.34 2.4237 0.005 **
## + Phosphorus 1 873.44 2.3340 0.005 **
## + Barium 1 873.62 2.1603 0.005 **
## + Nickel 1 873.63 2.1424 0.005 **
## + Moisture_Content 1 873.65 2.1231 0.005 **
## + Plot 1 873.71 2.0722 0.005 **
## + Manganese 1 873.73 2.0526 0.005 **
## + Nitrate_Nitrite 1 873.77 2.0084 0.005 **
## + Zinc 1 873.94 1.8456 0.005 **
## + Ammonia 1 873.95 1.8361 0.005 **
## + Calcium 1 874.00 1.7954 0.005 **
## + Iron 1 874.06 1.7324 0.005 **
## + Sulfur 1 874.08 1.7112 0.005 **
## + Copper 1 874.13 1.6675 0.005 **
## + grain_yield 1 874.17 1.6277 0.005 **
## + Cobalt 1 873.91 1.8798 0.015 *
## + pH 1 874.58 1.2350 0.130
## + Pea_variety 6 879.71 0.9109 0.850
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Tillage + prev_crop
##
## Df AIC F Pr(>F)
## - prev_crop 3 873.90 2.5825 0.005 **
## - Tillage 2 875.14 3.4833 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + irrgation 1 870.32 3.1471 0.005 **
## + Organic_Matter 1 870.67 2.8207 0.005 **
## + Ammonia 1 870.86 2.6513 0.005 **
## + total_precip_irr 1 871.34 2.2116 0.005 **
## + grain_yield 1 871.34 2.2108 0.005 **
## + Magnesium 1 871.34 2.2086 0.005 **
## + Barium 1 871.51 2.0507 0.005 **
## + season_precip 1 871.56 2.0064 0.005 **
## + Nickel 1 871.58 1.9892 0.005 **
## + Manganese 1 871.63 1.9468 0.005 **
## + pH 1 871.84 1.7541 0.005 **
## + Nitrate_Nitrite 1 871.63 1.9486 0.010 **
## + Cobalt 1 871.63 1.9418 0.010 **
## + Av_Potassium 1 871.94 1.6614 0.010 **
## + Calcium 1 872.01 1.6001 0.010 **
## + Zinc 1 871.95 1.6604 0.015 *
## + Copper 1 872.09 1.5325 0.020 *
## + Iron 1 871.81 1.7800 0.030 *
## + Av_Phosphorus 1 872.28 1.3598 0.045 *
## + Phosphorus 1 872.28 1.3633 0.120
## + Moisture_Content 1 872.41 1.2435 0.155
## + Sulfur 1 872.55 1.1177 0.340
## + Pea_variety 6 876.76 0.9773 0.570
## + Plot 0 871.82
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Tillage + prev_crop + irrgation
##
## Df AIC F Pr(>F)
## - irrgation 1 871.82 3.1471 0.005 **
## - prev_crop 3 872.87 2.6883 0.005 **
## - Tillage 2 872.96 3.0756 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Ammonia 1 870.03 1.9949 0.005 **
## + Manganese 1 870.12 1.9138 0.005 **
## + season_precip 1 870.13 1.9074 0.005 **
## + total_precip_irr 1 870.13 1.9074 0.005 **
## + grain_yield 1 870.20 1.8412 0.005 **
## + pH 1 870.27 1.7801 0.005 **
## + Av_Potassium 1 870.41 1.6562 0.005 **
## + Zinc 1 870.45 1.6182 0.005 **
## + Nitrate_Nitrite 1 870.01 2.0078 0.010 **
## + Organic_Matter 1 870.41 1.6546 0.010 **
## + Nickel 1 870.57 1.5191 0.010 **
## + Magnesium 1 870.54 1.5390 0.030 *
## + Copper 1 870.59 1.5016 0.040 *
## + Av_Phosphorus 1 870.70 1.3994 0.060 .
## + Cobalt 1 870.40 1.6632 0.090 .
## + Iron 1 870.65 1.4446 0.105
## + Moisture_Content 1 870.76 1.3526 0.115
## + Barium 1 870.95 1.1863 0.165
## + Phosphorus 1 870.83 1.2836 0.220
## + Sulfur 1 871.08 1.0655 0.375
## + Pea_variety 6 874.92 1.0036 0.520
## + Calcium 1 871.27 0.9011 0.685
## + Plot 0 870.32
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Tillage + prev_crop + irrgation + Ammonia
##
## Df AIC F Pr(>F)
## - Ammonia 1 870.32 1.9949 0.005 **
## - irrgation 1 870.86 2.4742 0.005 **
## - Tillage 2 872.15 2.7622 0.005 **
## - prev_crop 3 872.66 2.6590 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + season_precip 1 869.61 2.0562 0.005 **
## + total_precip_irr 1 869.61 2.0562 0.005 **
## + Nitrate_Nitrite 1 869.62 2.0545 0.005 **
## + pH 1 869.82 1.8745 0.005 **
## + Manganese 1 869.88 1.8270 0.005 **
## + Av_Potassium 1 870.01 1.7141 0.010 **
## + Organic_Matter 1 870.04 1.6910 0.010 **
## + Zinc 1 870.10 1.6350 0.010 **
## + grain_yield 1 870.11 1.6301 0.010 **
## + Nickel 1 870.13 1.6135 0.020 *
## + Copper 1 870.20 1.5505 0.025 *
## + Magnesium 1 870.33 1.4337 0.055 .
## + Iron 1 870.39 1.3872 0.135
## + Cobalt 1 870.43 1.3555 0.220
## + Av_Phosphorus 1 870.75 1.0745 0.250
## + Pea_variety 6 874.34 1.0205 0.365
## + Phosphorus 1 870.85 0.9908 0.435
## + Barium 1 870.85 0.9879 0.510
## + Sulfur 1 870.95 0.9063 0.625
## + Calcium 1 870.93 0.9229 0.655
## + Moisture_Content 1 871.02 0.8480 0.730
## + Plot 0 870.03
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Tillage + prev_crop + irrgation + Ammonia + season_precip
##
## Df AIC F Pr(>F)
## - season_precip 1 870.03 2.0562 0.005 **
## - Ammonia 1 870.13 2.1422 0.005 **
## - irrgation 1 870.58 2.5391 0.005 **
## - Tillage 2 871.49 2.5867 0.005 **
## - prev_crop 3 871.65 2.4051 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + grain_yield 1 869.54 1.7222 0.010 **
## + Nitrate_Nitrite 1 868.86 2.3031 0.020 *
## + Cobalt 1 869.93 1.3919 0.160
## + Manganese 1 869.86 1.4507 0.170
## + Av_Phosphorus 1 870.26 1.1212 0.270
## + Pea_variety 6 873.52 1.0515 0.310
## + Magnesium 1 870.30 1.0830 0.355
## + Iron 1 870.31 1.0795 0.360
## + pH 1 870.37 1.0279 0.415
## + Av_Potassium 1 870.39 1.0079 0.460
## + Phosphorus 1 870.39 1.0080 0.470
## + Barium 1 870.41 0.9929 0.515
## + Zinc 1 870.53 0.8910 0.575
## + Organic_Matter 1 870.50 0.9177 0.610
## + Nickel 1 870.60 0.8379 0.645
## + Sulfur 1 870.54 0.8863 0.660
## + Copper 1 870.61 0.8285 0.680
## + Moisture_Content 1 870.68 0.7700 0.800
## + Calcium 1 870.82 0.6544 0.915
## + Plot 0 869.61
## + total_precip_irr 0 869.61
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Tillage + prev_crop + irrgation + Ammonia + season_precip + grain_yield
##
## Df AIC F Pr(>F)
## - Ammonia 1 869.37 1.5161 0.015 *
## - grain_yield 1 869.61 1.7222 0.005 **
## - season_precip 1 870.11 2.1400 0.005 **
## - irrgation 1 870.60 2.5682 0.005 **
## - Tillage 2 871.60 2.6131 0.005 **
## - prev_crop 3 871.71 2.3952 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Nitrate_Nitrite 1 868.62 2.3907 0.010 **
## + Av_Phosphorus 1 870.05 1.2029 0.160
## + Pea_variety 6 873.33 1.0408 0.290
## + Cobalt 1 870.30 0.9985 0.415
## + Av_Potassium 1 870.34 0.9699 0.420
## + Iron 1 870.36 0.9536 0.425
## + Manganese 1 870.25 1.0415 0.465
## + pH 1 870.30 1.0004 0.505
## + Phosphorus 1 870.35 0.9605 0.535
## + Zinc 1 870.44 0.8877 0.540
## + Magnesium 1 870.46 0.8682 0.605
## + Barium 1 870.42 0.9015 0.620
## + Copper 1 870.49 0.8450 0.685
## + Nickel 1 870.49 0.8494 0.740
## + Sulfur 1 870.62 0.7421 0.830
## + Organic_Matter 1 870.62 0.7391 0.850
## + Moisture_Content 1 870.65 0.7188 0.890
## + Calcium 1 870.70 0.6781 0.905
## + Plot 0 869.54
## + total_precip_irr 0 869.54
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Tillage + prev_crop + irrgation + Ammonia + season_precip + grain_yield + Nitrate_Nitrite
##
## Df AIC F Pr(>F)
## - Ammonia 1 868.61 1.6173 0.015 *
## - grain_yield 1 868.86 1.8213 0.005 **
## - season_precip 1 869.33 2.2157 0.005 **
## - Nitrate_Nitrite 1 869.54 2.3907 0.005 **
## - irrgation 1 869.85 2.6529 0.005 **
## - Tillage 2 870.93 2.6640 0.005 **
## - prev_crop 3 871.16 2.4551 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Av_Phosphorus 1 869.14 1.1639 0.225
## + Pea_variety 6 872.23 1.0371 0.365
## + Cobalt 1 869.31 1.0309 0.370
## + Manganese 1 869.27 1.0642 0.380
## + Av_Potassium 1 869.31 1.0299 0.390
## + pH 1 869.33 1.0128 0.435
## + Phosphorus 1 869.31 1.0293 0.450
## + Iron 1 869.37 0.9871 0.485
## + Zinc 1 869.47 0.9046 0.525
## + Magnesium 1 869.45 0.9228 0.580
## + Barium 1 869.46 0.9115 0.595
## + Copper 1 869.46 0.9110 0.605
## + Nickel 1 869.55 0.8407 0.780
## + Moisture_Content 1 869.71 0.7150 0.880
## + Calcium 1 869.73 0.6989 0.880
## + Sulfur 1 869.67 0.7431 0.890
## + Organic_Matter 1 869.68 0.7370 0.905
## + Plot 0 868.62
## + total_precip_irr 0 868.62m1_nifH_cca_chem <- cca(abund_tablenifH ~ Organic_Matter + Moisture_Content + Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium + pH + Barium + Calcium + Cobalt + Copper + Iron + Magnesium + Manganese + Nickel + Phosphorus + Sulfur + Zinc , meta3)
m0_nifH_cca_chem <- cca(abund_tablenifH ~ 1, meta3)
m1_nifH_cca_chem## Call: cca(formula = abund_tablenifH ~ Organic_Matter +
## Moisture_Content + Nitrate_Nitrite + Ammonia + Av_Phosphorus +
## Av_Potassium + pH + Barium + Calcium + Cobalt + Copper + Iron +
## Magnesium + Manganese + Nickel + Phosphorus + Sulfur + Zinc, data
## = meta3)
##
## Inertia Proportion Rank
## Total 10.7659 1.0000
## Constrained 5.2066 0.4836 18
## Unconstrained 5.5593 0.5164 35
## Inertia is scaled Chi-square
##
## Eigenvalues for constrained axes:
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7 CCA8 CCA9 CCA10
## 0.7916 0.6368 0.5430 0.5237 0.4736 0.4137 0.3403 0.2283 0.2008 0.1907
## CCA11 CCA12 CCA13 CCA14 CCA15 CCA16 CCA17 CCA18
## 0.1609 0.1430 0.1267 0.1085 0.0994 0.0838 0.0794 0.0625
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.4573 0.4358 0.3745 0.3252 0.2860 0.2803 0.2665 0.2449
## (Showed only 8 of all 35 unconstrained eigenvalues)m0_nifH_cca_chem## Call: cca(formula = abund_tablenifH ~ 1, data = meta3)
##
## Inertia Rank
## Total 10.77
## Unconstrained 10.77 53
## Inertia is scaled Chi-square
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.8189 0.6921 0.6401 0.5860 0.5584 0.5037 0.4818 0.4700
## (Showed only 8 of all 53 unconstrained eigenvalues)model_nifH_cca_chem <-ordistep(m0_nifH_cca_chem, scope=formula(m1_nifH_cca_chem))##
## Start: abund_tablenifH ~ 1
##
## Df AIC F Pr(>F)
## + Organic_Matter 1 874.75 3.5957 0.005 **
## + Nitrate_Nitrite 1 875.40 2.9283 0.005 **
## + Zinc 1 875.68 2.6463 0.005 **
## + Av_Potassium 1 875.71 2.6178 0.005 **
## + Nickel 1 875.88 2.4512 0.005 **
## + Cobalt 1 875.92 2.4074 0.005 **
## + Barium 1 875.93 2.3929 0.005 **
## + Iron 1 875.93 2.3923 0.005 **
## + Copper 1 875.99 2.3375 0.005 **
## + Magnesium 1 876.00 2.3264 0.005 **
## + Av_Phosphorus 1 876.01 2.3202 0.005 **
## + Phosphorus 1 876.07 2.2549 0.005 **
## + Moisture_Content 1 876.18 2.1464 0.005 **
## + Manganese 1 876.18 2.1462 0.005 **
## + Ammonia 1 876.47 1.8528 0.005 **
## + Calcium 1 876.57 1.7584 0.005 **
## + Sulfur 1 876.47 1.8514 0.010 **
## + pH 1 876.84 1.4897 0.030 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter
##
## Df AIC F Pr(>F)
## - Organic_Matter 1 876.36 3.5957 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Nitrate_Nitrite 1 873.56 3.1042 0.005 **
## + Zinc 1 873.97 2.6966 0.005 **
## + Barium 1 874.16 2.5125 0.005 **
## + Magnesium 1 874.20 2.4698 0.005 **
## + Cobalt 1 874.21 2.4614 0.005 **
## + Copper 1 874.22 2.4475 0.005 **
## + Av_Potassium 1 874.23 2.4431 0.005 **
## + Iron 1 874.24 2.4244 0.005 **
## + Av_Phosphorus 1 874.27 2.3965 0.005 **
## + Phosphorus 1 874.34 2.3325 0.005 **
## + Nickel 1 874.36 2.3098 0.005 **
## + Manganese 1 874.39 2.2789 0.005 **
## + Moisture_Content 1 874.48 2.1872 0.005 **
## + Sulfur 1 874.63 2.0475 0.005 **
## + Calcium 1 874.79 1.8862 0.005 **
## + Ammonia 1 874.88 1.7998 0.005 **
## + pH 1 875.14 1.5430 0.010 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite
##
## Df AIC F Pr(>F)
## - Nitrate_Nitrite 1 874.75 3.1042 0.005 **
## - Organic_Matter 1 875.40 3.7617 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Av_Potassium 1 872.79 2.6309 0.005 **
## + Magnesium 1 872.92 2.5090 0.005 **
## + Av_Phosphorus 1 872.98 2.4507 0.005 **
## + Phosphorus 1 872.99 2.4422 0.005 **
## + Moisture_Content 1 873.15 2.2870 0.005 **
## + Manganese 1 873.22 2.2192 0.005 **
## + Zinc 1 873.36 2.0786 0.005 **
## + Sulfur 1 873.39 2.0553 0.005 **
## + Copper 1 873.50 1.9447 0.005 **
## + Cobalt 1 873.54 1.9069 0.005 **
## + Calcium 1 873.55 1.8973 0.005 **
## + Barium 1 873.61 1.8369 0.005 **
## + Ammonia 1 873.65 1.8007 0.005 **
## + Iron 1 873.67 1.7782 0.005 **
## + Nickel 1 873.84 1.6164 0.005 **
## + pH 1 873.98 1.4887 0.020 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium
##
## Df AIC F Pr(>F)
## - Av_Potassium 1 873.56 2.6309 0.005 **
## - Nitrate_Nitrite 1 874.23 3.2820 0.005 **
## - Organic_Matter 1 874.51 3.5682 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Av_Phosphorus 1 872.25 2.3577 0.005 **
## + Magnesium 1 872.36 2.2606 0.005 **
## + Manganese 1 872.56 2.0639 0.005 **
## + Sulfur 1 872.61 2.0249 0.005 **
## + Calcium 1 872.63 2.0019 0.005 **
## + Ammonia 1 872.75 1.8850 0.005 **
## + Barium 1 872.80 1.8420 0.005 **
## + Nickel 1 872.90 1.7454 0.005 **
## + pH 1 872.93 1.7153 0.005 **
## + Copper 1 872.96 1.6867 0.005 **
## + Zinc 1 872.98 1.6722 0.005 **
## + Phosphorus 1 873.06 1.5945 0.005 **
## + Cobalt 1 872.89 1.7581 0.010 **
## + Moisture_Content 1 873.18 1.4875 0.015 *
## + Iron 1 873.16 1.4996 0.035 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium + Av_Phosphorus
##
## Df AIC F Pr(>F)
## - Av_Phosphorus 1 872.79 2.3577 0.005 **
## - Av_Potassium 1 872.98 2.5341 0.005 **
## - Nitrate_Nitrite 1 873.69 3.2208 0.005 **
## - Organic_Matter 1 874.15 3.6662 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Manganese 1 871.90 2.1402 0.005 **
## + Magnesium 1 872.03 2.0192 0.005 **
## + Sulfur 1 872.11 1.9489 0.005 **
## + Calcium 1 872.19 1.8727 0.005 **
## + Barium 1 872.29 1.7817 0.005 **
## + Ammonia 1 872.36 1.7145 0.005 **
## + Zinc 1 872.38 1.6967 0.005 **
## + Copper 1 872.44 1.6399 0.005 **
## + Phosphorus 1 872.48 1.6044 0.005 **
## + Cobalt 1 872.28 1.7883 0.010 **
## + Nickel 1 872.47 1.6161 0.010 **
## + Iron 1 872.54 1.5510 0.010 **
## + pH 1 872.58 1.5087 0.025 *
## + Moisture_Content 1 872.82 1.2883 0.080 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium + Av_Phosphorus + Manganese
##
## Df AIC F Pr(>F)
## - Manganese 1 872.25 2.1402 0.005 **
## - Av_Phosphorus 1 872.56 2.4287 0.005 **
## - Av_Potassium 1 872.74 2.5890 0.005 **
## - Nitrate_Nitrite 1 873.12 2.9469 0.005 **
## - Organic_Matter 1 873.85 3.6433 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Barium 1 871.44 2.1898 0.005 **
## + Sulfur 1 871.71 1.9435 0.005 **
## + Magnesium 1 871.76 1.9029 0.005 **
## + Zinc 1 871.79 1.8710 0.005 **
## + Cobalt 1 871.82 1.8413 0.005 **
## + Iron 1 871.91 1.7662 0.005 **
## + Copper 1 871.91 1.7627 0.005 **
## + Calcium 1 871.92 1.7571 0.005 **
## + Nickel 1 872.03 1.6558 0.005 **
## + Phosphorus 1 872.05 1.6338 0.005 **
## + pH 1 872.11 1.5817 0.005 **
## + Ammonia 1 872.20 1.5044 0.010 **
## + Moisture_Content 1 872.38 1.3407 0.050 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium + Av_Phosphorus + Manganese + Barium
##
## Df AIC F Pr(>F)
## - Barium 1 871.90 2.1898 0.005 **
## - Av_Phosphorus 1 872.08 2.3589 0.005 **
## - Nitrate_Nitrite 1 872.21 2.4735 0.005 **
## - Av_Potassium 1 872.29 2.5438 0.005 **
## - Manganese 1 872.29 2.5440 0.005 **
## - Organic_Matter 1 873.10 3.2928 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Cobalt 1 871.21 1.9353 0.005 **
## + Copper 1 871.30 1.8613 0.005 **
## + Sulfur 1 871.30 1.8589 0.005 **
## + Iron 1 871.31 1.8503 0.005 **
## + Magnesium 1 871.35 1.8125 0.005 **
## + Zinc 1 871.42 1.7495 0.005 **
## + Calcium 1 871.43 1.7405 0.005 **
## + Phosphorus 1 871.56 1.6326 0.005 **
## + Nickel 1 871.65 1.5524 0.005 **
## + pH 1 871.56 1.6330 0.010 **
## + Moisture_Content 1 871.84 1.3860 0.015 *
## + Ammonia 1 871.84 1.3831 0.025 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium + Av_Phosphorus + Manganese + Barium + Cobalt
##
## Df AIC F Pr(>F)
## - Cobalt 1 871.44 1.9353 0.005 **
## - Barium 1 871.82 2.2773 0.005 **
## - Av_Phosphorus 1 871.87 2.3225 0.005 **
## - Nitrate_Nitrite 1 871.91 2.3545 0.005 **
## - Manganese 1 871.92 2.3678 0.005 **
## - Av_Potassium 1 872.23 2.6409 0.005 **
## - Organic_Matter 1 873.08 3.4110 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Sulfur 1 870.86 2.0053 0.005 **
## + Calcium 1 870.96 1.9217 0.005 **
## + Nickel 1 871.19 1.7216 0.005 **
## + pH 1 871.48 1.4699 0.010 **
## + Phosphorus 1 871.57 1.3904 0.025 *
## + Moisture_Content 1 871.58 1.3814 0.025 *
## + Ammonia 1 871.61 1.3613 0.040 *
## + Zinc 1 871.68 1.2995 0.060 .
## + Magnesium 1 871.71 1.2731 0.070 .
## + Copper 1 871.87 1.1341 0.215
## + Iron 1 872.01 1.0150 0.530
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium + Av_Phosphorus + Manganese + Barium + Cobalt + Sulfur
##
## Df AIC F Pr(>F)
## - Sulfur 1 871.21 2.0053 0.005 **
## - Cobalt 1 871.30 2.0804 0.005 **
## - Barium 1 871.39 2.1625 0.005 **
## - Av_Phosphorus 1 871.44 2.1990 0.005 **
## - Manganese 1 871.57 2.3147 0.005 **
## - Nitrate_Nitrite 1 871.67 2.4073 0.005 **
## - Av_Potassium 1 872.03 2.7182 0.005 **
## - Organic_Matter 1 872.99 3.5736 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Magnesium 1 871.12 1.4435 0.010 **
## + Nickel 1 871.12 1.4399 0.040 *
## + pH 1 871.36 1.2417 0.090 .
## + Copper 1 871.38 1.2255 0.125
## + Phosphorus 1 871.35 1.2508 0.130
## + Zinc 1 871.41 1.2013 0.145
## + Ammonia 1 871.53 1.0985 0.305
## + Moisture_Content 1 871.63 1.0163 0.485
## + Iron 1 871.66 0.9862 0.520
## + Calcium 1 871.93 0.7671 0.835
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium + Av_Phosphorus + Manganese + Barium + Cobalt + Sulfur + Magnesium
##
## Df AIC F Pr(>F)
## - Cobalt 1 870.55 1.1831 0.225
## - Magnesium 1 870.86 1.4435 0.030 *
## - Av_Phosphorus 1 870.94 1.5081 0.005 **
## - Barium 1 871.20 1.7285 0.005 **
## - Sulfur 1 871.71 2.1626 0.005 **
## - Manganese 1 871.72 2.1755 0.005 **
## - Nitrate_Nitrite 1 871.80 2.2421 0.005 **
## - Organic_Matter 1 871.94 2.3606 0.005 **
## - Av_Potassium 1 872.38 2.7380 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium + Av_Phosphorus + Manganese + Barium + Sulfur + Magnesium
##
## Df AIC F Pr(>F)
## + Phosphorus 1 870.39 1.7956 0.005 **
## + Iron 1 870.99 1.2859 0.155
## + Ammonia 1 871.17 1.1384 0.210
## + pH 1 871.19 1.1215 0.230
## + Cobalt 1 871.12 1.1831 0.245
## + Nickel 1 871.20 1.1160 0.310
## + Zinc 1 871.44 0.9148 0.650
## + Moisture_Content 1 871.48 0.8821 0.675
## + Copper 1 871.45 0.9088 0.710
## + Calcium 1 871.50 0.8643 0.740
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium + Av_Phosphorus + Manganese + Barium + Sulfur + Magnesium + Phosphorus
##
## Df AIC F Pr(>F)
## - Av_Phosphorus 1 870.15 1.4610 0.025 *
## - Av_Potassium 1 870.27 1.5617 0.010 **
## - Phosphorus 1 870.55 1.7956 0.005 **
## - Magnesium 1 870.72 1.9382 0.005 **
## - Nitrate_Nitrite 1 870.93 2.1213 0.005 **
## - Sulfur 1 871.03 2.2060 0.005 **
## - Manganese 1 871.07 2.2387 0.005 **
## - Barium 1 871.10 2.2666 0.005 **
## - Organic_Matter 1 871.48 2.5897 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + pH 1 870.92 1.1871 0.180
## + Ammonia 1 870.94 1.1697 0.215
## + Nickel 1 871.03 1.0927 0.275
## + Cobalt 1 871.04 1.0916 0.275
## + Iron 1 871.03 1.0935 0.315
## + Zinc 1 871.13 1.0140 0.445
## + Calcium 1 871.23 0.9362 0.510
## + Moisture_Content 1 871.28 0.8899 0.735
## + Copper 1 871.33 0.8507 0.775model_nifH_cca_chem$anova## Df AIC F Pr(>F)
## + Organic_Matter 1 874.75 3.5957 0.005 **
## + Nitrate_Nitrite 1 873.56 3.1042 0.005 **
## + Av_Potassium 1 872.79 2.6309 0.005 **
## + Av_Phosphorus 1 872.25 2.3577 0.005 **
## + Manganese 1 871.90 2.1402 0.005 **
## + Barium 1 871.44 2.1898 0.005 **
## + Cobalt 1 871.21 1.9353 0.005 **
## + Sulfur 1 870.86 2.0053 0.005 **
## + Magnesium 1 871.12 1.4435 0.010 **
## - Cobalt 1 870.55 1.1831 0.225
## + Phosphorus 1 870.39 1.7956 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1RDA is a linear cca
m1_nifH_rda <- rda(abund_tablenifH ~ ., meta3)
m0_nifH_rda <- rda(abund_tablenifH ~ 1, meta3)
m1_nifH_rda## Call: rda(formula = abund_tablenifH ~ Site + Pea_variety + Plot +
## season_precip + irrgation + total_precip_irr + Tillage + prev_crop
## + grain_yield + Organic_Matter + Moisture_Content +
## Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium + pH +
## Barium + Calcium + Cobalt + Copper + Iron + Magnesium + Manganese
## + Nickel + Phosphorus + Sulfur + Zinc, data = meta3)
##
## Inertia Proportion Rank
## Total 1.258e+11 1.000e+00
## Constrained 9.747e+10 7.745e-01 33
## Unconstrained 2.837e+10 2.255e-01 20
## Inertia is variance
## Some constraints were aliased because they were collinear (redundant)
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## 16353665417 15356453441 12209872418 10446884520 7322044854 5776558730
## RDA7 RDA8 RDA9 RDA10 RDA11 RDA12
## 4833284412 4476474731 3654563170 3048908811 2374942944 1890413915
## RDA13 RDA14 RDA15 RDA16 RDA17 RDA18
## 1488069677 1072841995 989712375 941967373 836809708 626843871
## RDA19 RDA20 RDA21 RDA22 RDA23 RDA24
## 573856607 508536865 403729752 376714333 325054973 282465322
## RDA25 RDA26 RDA27 RDA28 RDA29 RDA30
## 238159099 207287067 185699351 163253845 127967489 116333192
## RDA31 RDA32 RDA33
## 95605085 84351951 75938358
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6
## 5599447125 4342938458 3880332013 3368259818 2027561149 1483067201
## PC7 PC8
## 1310478544 1156592863
## (Showed only 8 of all 20 unconstrained eigenvalues)m0_nifH_rda## Call: rda(formula = abund_tablenifH ~ 1, data = meta3)
##
## Inertia Rank
## Total 1.258e+11
## Unconstrained 1.258e+11 53
## Inertia is variance
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6
## 18388165460 17017988044 13929669751 13228540294 10298920752 7474983345
## PC7 PC8
## 6353906098 5386444205
## (Showed only 8 of all 53 unconstrained eigenvalues)plot(m1_nifH_rda)
model_rda_nifH <-ordiR2step(m0_nifH_rda, scope=formula(m1_nifH_rda))## Step: R2.adj= 0
## Call: abund_tablenifH ~ 1
##
## R2.adjusted
## <All variables> 0.40253072
## + Site 0.39284983
## + prev_crop 0.14946826
## + Tillage 0.09118753
## + total_precip_irr 0.06245835
## + Av_Phosphorus 0.06052058
## + Phosphorus 0.05195778
## + Av_Potassium 0.05112278
## + Magnesium 0.04894792
## + irrgation 0.04678398
## + Plot 0.04077308
## + Moisture_Content 0.04047656
## + season_precip 0.03754349
## + Barium 0.03466802
## + Manganese 0.03434338
## + Zinc 0.03212237
## + Organic_Matter 0.03177594
## + Ammonia 0.02812448
## + Nitrate_Nitrite 0.02759831
## + Nickel 0.02460517
## + Copper 0.02201680
## + Sulfur 0.02175055
## + Calcium 0.01980531
## + grain_yield 0.01948985
## + Iron 0.01591920
## + pH 0.01585647
## + Cobalt 0.01268035
## <none> 0.00000000
## + Pea_variety -0.02472576
##
## Df AIC F Pr(>F)
## + Site 8 1361.4 5.2866 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.3928498
## Call: abund_tablenifH ~ Site
##
## R2.adjusted
## + Pea_variety 0.4053417
## <All variables> 0.4025307
## + Ammonia 0.3971916
## + Nitrate_Nitrite 0.3945793
## <none> 0.3928498
## + Plot 0.3928498
## + season_precip 0.3928498
## + irrgation 0.3928498
## + total_precip_irr 0.3928498
## + Tillage 0.3928498
## + prev_crop 0.3928498
## + Av_Potassium 0.3920019
## + Organic_Matter 0.3905573
## + Av_Phosphorus 0.3888507
## + Sulfur 0.3888353
## + pH 0.3883561
## + Calcium 0.3881788
## + Moisture_Content 0.3878418
## + Manganese 0.3866273
## + grain_yield 0.3860617
## + Magnesium 0.3853474
## + Nickel 0.3852829
## + Cobalt 0.3846202
## + Copper 0.3845997
## + Iron 0.3842913
## + Barium 0.3842156
## + Phosphorus 0.3841708
## + Zinc 0.3838492Chemistry only model
m1_nifH_rda_chem <- rda(abund_tablenifH ~ Organic_Matter + Moisture_Content + Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium + pH + Barium + Calcium + Cobalt + Copper + Iron + Magnesium + Manganese + Nickel + Phosphorus + Sulfur + Zinc , meta3)
m0_nifH_rda_chem <- rda(abund_tablenifH ~ 1, meta3)
m1_nifH_rda_chem## Call: rda(formula = abund_tablenifH ~ Organic_Matter +
## Moisture_Content + Nitrate_Nitrite + Ammonia + Av_Phosphorus +
## Av_Potassium + pH + Barium + Calcium + Cobalt + Copper + Iron +
## Magnesium + Manganese + Nickel + Phosphorus + Sulfur + Zinc, data
## = meta3)
##
## Inertia Proportion Rank
## Total 1.258e+11 1.000e+00
## Constrained 6.893e+10 5.478e-01 18
## Unconstrained 5.690e+10 4.522e-01 35
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## 14804112412 13108358521 10437172065 8136887370 5686929136 3911450523
## RDA7 RDA8 RDA9 RDA10 RDA11 RDA12
## 3487763724 2618476736 1896564759 1302627898 875978909 653782995
## RDA13 RDA14 RDA15 RDA16 RDA17 RDA18
## 566194235 425743412 312506701 290536181 223053545 193225475
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6
## 11018590525 7155476437 5939185648 5422891779 4017288576 3016129382
## PC7 PC8
## 2886770730 2071617511
## (Showed only 8 of all 35 unconstrained eigenvalues)m0_nifH_rda_chem## Call: rda(formula = abund_tablenifH ~ 1, data = meta3)
##
## Inertia Rank
## Total 1.258e+11
## Unconstrained 1.258e+11 53
## Inertia is variance
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6
## 18388165460 17017988044 13929669751 13228540294 10298920752 7474983345
## PC7 PC8
## 6353906098 5386444205
## (Showed only 8 of all 53 unconstrained eigenvalues)plot(m1_nifH_rda_chem)
model_rda_chem_nifH <-ordiR2step(m0_nifH_rda_chem, scope=formula(m1_nifH_rda_chem))## Step: R2.adj= 0
## Call: abund_tablenifH ~ 1
##
## R2.adjusted
## <All variables> 0.31521855
## + Av_Phosphorus 0.06052058
## + Phosphorus 0.05195778
## + Av_Potassium 0.05112278
## + Magnesium 0.04894792
## + Moisture_Content 0.04047656
## + Barium 0.03466802
## + Manganese 0.03434338
## + Zinc 0.03212237
## + Organic_Matter 0.03177594
## + Ammonia 0.02812448
## + Nitrate_Nitrite 0.02759831
## + Nickel 0.02460517
## + Copper 0.02201680
## + Sulfur 0.02175055
## + Calcium 0.01980531
## + Iron 0.01591920
## + pH 0.01585647
## + Cobalt 0.01268035
## <none> 0.00000000
##
## Df AIC F Pr(>F)
## + Av_Phosphorus 1 1378.7 4.4142 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.06052058
## Call: abund_tablenifH ~ Av_Phosphorus
##
## R2.adjusted
## <All variables> 0.31521855
## + Av_Potassium 0.10253413
## + Phosphorus 0.09795318
## + Zinc 0.09505177
## + Manganese 0.09455271
## + Organic_Matter 0.09324882
## + Ammonia 0.09270795
## + Nickel 0.08781964
## + Magnesium 0.08490752
## + Barium 0.08424954
## + Moisture_Content 0.08392402
## + Nitrate_Nitrite 0.08340692
## + Sulfur 0.08244671
## + Copper 0.08171114
## + Iron 0.08134433
## + Calcium 0.07884137
## + Cobalt 0.07764768
## + pH 0.07328876
## <none> 0.06052058
##
## Df AIC F Pr(>F)
## + Av_Potassium 1 1377.2 3.4343 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.1025341
## Call: abund_tablenifH ~ Av_Phosphorus + Av_Potassium
##
## R2.adjusted
## <All variables> 0.3152185
## + Organic_Matter 0.1426010
## + Zinc 0.1346765
## + Manganese 0.1346442
## + Sulfur 0.1321684
## + Nickel 0.1314652
## + Phosphorus 0.1297865
## + Calcium 0.1276961
## + Magnesium 0.1269648
## + Nitrate_Nitrite 0.1268511
## + Ammonia 0.1264398
## + Barium 0.1263189
## + Copper 0.1232790
## + Moisture_Content 0.1227996
## + pH 0.1220798
## + Iron 0.1192709
## + Cobalt 0.1140726
## <none> 0.1025341
##
## Df AIC F Pr(>F)
## + Organic_Matter 1 1375.7 3.3833 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.142601
## Call: abund_tablenifH ~ Av_Phosphorus + Av_Potassium + Organic_Matter
##
## R2.adjusted
## <All variables> 0.3152185
## + Barium 0.1806009
## + Manganese 0.1778531
## + Zinc 0.1739765
## + Ammonia 0.1728053
## + Sulfur 0.1698627
## + Calcium 0.1691275
## + Magnesium 0.1685205
## + Nitrate_Nitrite 0.1637934
## + Nickel 0.1637017
## + Copper 0.1636724
## + pH 0.1635979
## + Phosphorus 0.1610160
## + Iron 0.1572052
## + Moisture_Content 0.1536428
## + Cobalt 0.1526631
## <none> 0.1426010
##
## Df AIC F Pr(>F)
## + Barium 1 1374.1 3.3188 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.1806009
## Call: abund_tablenifH ~ Av_Phosphorus + Av_Potassium + Organic_Matter + Barium
##
## R2.adjusted
## <All variables> 0.3152185
## + Manganese 0.2425367
## + Sulfur 0.2152786
## + Calcium 0.2146288
## + Magnesium 0.2075492
## + Cobalt 0.2059481
## + Iron 0.2045766
## + Zinc 0.2041698
## + Nickel 0.2040281
## + pH 0.2023023
## + Ammonia 0.2016233
## + Copper 0.2015660
## + Phosphorus 0.2011897
## + Moisture_Content 0.1920917
## + Nitrate_Nitrite 0.1880401
## <none> 0.1806009
##
## Df AIC F Pr(>F)
## + Manganese 1 1370.8 5.0066 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.2425367
## Call: abund_tablenifH ~ Av_Phosphorus + Av_Potassium + Organic_Matter + Barium + Manganese
##
## R2.adjusted
## <All variables> 0.3152185
## + Iron 0.2741671
## + Magnesium 0.2715470
## + Zinc 0.2703388
## + Sulfur 0.2688961
## + Copper 0.2683106
## + Cobalt 0.2675120
## + pH 0.2661270
## + Calcium 0.2650136
## + Phosphorus 0.2645607
## + Nickel 0.2626842
## + Moisture_Content 0.2555896
## + Ammonia 0.2539120
## + Nitrate_Nitrite 0.2511329
## <none> 0.2425367
##
## Df AIC F Pr(>F)
## + Iron 1 1369.3 3.0917 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.2741671
## Call: abund_tablenifH ~ Av_Phosphorus + Av_Potassium + Organic_Matter + Barium + Manganese + Iron
##
## R2.adjusted
## <All variables> 0.3152185
## + Sulfur 0.3077588
## + Calcium 0.3034844
## + pH 0.2895106
## + Phosphorus 0.2873660
## + Moisture_Content 0.2869075
## + Ammonia 0.2867278
## + Nickel 0.2850464
## + Zinc 0.2833410
## + Magnesium 0.2802296
## + Nitrate_Nitrite 0.2787387
## + Cobalt 0.2743731
## <none> 0.2741671
## + Copper 0.2741048
##
## Df AIC F Pr(>F)
## + Sulfur 1 1367.6 3.2807 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.3077588
## Call: abund_tablenifH ~ Av_Phosphorus + Av_Potassium + Organic_Matter + Barium + Manganese + Iron + Sulfur
##
## R2.adjusted
## + Magnesium 0.3193868
## + Zinc 0.3191464
## <All variables> 0.3152185
## + Ammonia 0.3141646
## + pH 0.3132794
## + Phosphorus 0.3131892
## + Nitrate_Nitrite 0.3129262
## + Copper 0.3093376
## + Moisture_Content 0.3090885
## + Nickel 0.3081999
## <none> 0.3077588
## + Cobalt 0.3049041
## + Calcium 0.3007448model_rda_chem_nifH$anova## R2.adj Df AIC F Pr(>F)
## + Av_Phosphorus 0.060521 1 1378.7 4.4142 0.002 **
## + Av_Potassium 0.102534 1 1377.2 3.4343 0.002 **
## + Organic_Matter 0.142601 1 1375.7 3.3833 0.002 **
## + Barium 0.180601 1 1374.1 3.3188 0.002 **
## + Manganese 0.242537 1 1370.8 5.0066 0.002 **
## + Iron 0.274167 1 1369.3 3.0917 0.002 **
## + Sulfur 0.307759 1 1367.6 3.2807 0.002 **
## <All variables> 0.315219
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Our RDA chem model can explain 30% of the variance in the bacterial community.
CAP Ordination model building
m1_nifH_cap_chem<- capscale(abund_tablenifH ~ Organic_Matter + Moisture_Content + Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium + pH + Barium + Calcium + Cobalt + Copper + Iron + Magnesium + Manganese + Nickel + Phosphorus + Sulfur + Zinc , data = meta3, distance = "bray")
m0_nifH_cap_chem<- capscale(abund_distnifH ~ 1, meta3)
plot(m1_nifH_cap_chem)
model_cap_chem_nifH <-ordiR2step(m0_nifH_cap_chem, scope=formula(m1_nifH_cap_chem))## Step: R2.adj= 0
## Call: abund_distnifH ~ 1
##
## R2.adjusted
## <All variables> 0.42412042
## + Organic_Matter 0.08664823
## + Av_Potassium 0.06424058
## + Av_Phosphorus 0.05197708
## + Zinc 0.04843832
## + Nitrate_Nitrite 0.04757351
## + Nickel 0.04703479
## + Barium 0.04533845
## + Iron 0.04135991
## + Magnesium 0.04111715
## + Cobalt 0.03937854
## + Phosphorus 0.03739460
## + Manganese 0.03679948
## + Copper 0.03582041
## + Ammonia 0.03314279
## + Moisture_Content 0.02637968
## + Sulfur 0.02145207
## + pH 0.02106463
## + Calcium 0.02096101
## <none> 0.00000000
##
## Df AIC F Pr(>F)
## + Organic_Matter 1 159.43 5.8777 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.08664823
## Call: abund_distnifH ~ Organic_Matter
##
## R2.adjusted
## <All variables> 0.42412042
## + Barium 0.14270722
## + Av_Potassium 0.14238286
## + Av_Phosphorus 0.13935754
## + Nitrate_Nitrite 0.13484559
## + Zinc 0.13050213
## + Magnesium 0.12895812
## + Phosphorus 0.12745386
## + Manganese 0.12724465
## + Copper 0.12206419
## + Iron 0.12132345
## + Cobalt 0.12126082
## + Ammonia 0.11844446
## + Nickel 0.11807834
## + Sulfur 0.11669887
## + Moisture_Content 0.11655745
## + Calcium 0.11156606
## + pH 0.10979577
## <none> 0.08664823
##
## Df AIC F Pr(>F)
## + Barium 1 157.08 4.2814 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.1427072
## Call: abund_distnifH ~ Organic_Matter + Barium
##
## R2.adjusted
## <All variables> 0.4241204
## + Av_Potassium 0.2008529
## + Av_Phosphorus 0.1900137
## + Manganese 0.1880846
## + Phosphorus 0.1853480
## + Magnesium 0.1842902
## + Cobalt 0.1825912
## + Copper 0.1820141
## + Zinc 0.1809508
## + Iron 0.1802010
## + Ammonia 0.1754906
## + Moisture_Content 0.1750410
## + Nickel 0.1717581
## + Sulfur 0.1707810
## + Calcium 0.1678310
## + pH 0.1675754
## + Nitrate_Nitrite 0.1669577
## <none> 0.1427072
##
## Df AIC F Pr(>F)
## + Av_Potassium 1 154.35 4.5718 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.2008529
## Call: abund_distnifH ~ Organic_Matter + Barium + Av_Potassium
##
## R2.adjusted
## <All variables> 0.4241204
## + Manganese 0.2571847
## + Av_Phosphorus 0.2465866
## + Cobalt 0.2426324
## + Copper 0.2399983
## + Magnesium 0.2390988
## + Iron 0.2386871
## + Ammonia 0.2363363
## + Zinc 0.2360117
## + Nickel 0.2320877
## + pH 0.2311367
## + Sulfur 0.2301444
## + Calcium 0.2287892
## + Phosphorus 0.2236127
## + Nitrate_Nitrite 0.2227016
## + Moisture_Content 0.2141867
## <none> 0.2008529
##
## Df AIC F Pr(>F)
## + Manganese 1 151.47 4.637 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.2571847
## Call: abund_distnifH ~ Organic_Matter + Barium + Av_Potassium + Manganese
##
## R2.adjusted
## <All variables> 0.4241204
## + Av_Phosphorus 0.3048552
## + Copper 0.2999316
## + Magnesium 0.2990454
## + Iron 0.2975378
## + Zinc 0.2954468
## + Cobalt 0.2944697
## + pH 0.2893062
## + Nickel 0.2881294
## + Ammonia 0.2861456
## + Sulfur 0.2846872
## + Calcium 0.2815357
## + Phosphorus 0.2812081
## + Nitrate_Nitrite 0.2764049
## + Moisture_Content 0.2720061
## <none> 0.2571847
##
## Df AIC F Pr(>F)
## + Av_Phosphorus 1 148.93 4.2071 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.3048552
## Call: abund_distnifH ~ Organic_Matter + Barium + Av_Potassium + Manganese + Av_Phosphorus
##
## R2.adjusted
## <All variables> 0.4241204
## + Iron 0.3478017
## + Cobalt 0.3450564
## + Copper 0.3442042
## + Magnesium 0.3418955
## + Zinc 0.3394628
## + pH 0.3331523
## + Sulfur 0.3315043
## + Phosphorus 0.3300620
## + Nickel 0.3283710
## + Calcium 0.3275764
## + Nitrate_Nitrite 0.3265021
## + Moisture_Content 0.3154696
## + Ammonia 0.3126530
## <none> 0.3048552
##
## Df AIC F Pr(>F)
## + Iron 1 146.52 4.0021 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.3478017
## Call: abund_distnifH ~ Organic_Matter + Barium + Av_Potassium + Manganese + Av_Phosphorus + Iron
##
## R2.adjusted
## <All variables> 0.4241204
## + Sulfur 0.3885315
## + Calcium 0.3839583
## + Nickel 0.3693593
## + pH 0.3689149
## + Nitrate_Nitrite 0.3651795
## + Zinc 0.3638021
## + Phosphorus 0.3632509
## + Moisture_Content 0.3622189
## + Magnesium 0.3594648
## + Ammonia 0.3562554
## + Copper 0.3550172
## + Cobalt 0.3480682
## <none> 0.3478017
##
## Df AIC F Pr(>F)
## + Sulfur 1 144.06 3.9597 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.3885315
## Call: abund_distnifH ~ Organic_Matter + Barium + Av_Potassium + Manganese + Av_Phosphorus + Iron + Sulfur
##
## R2.adjusted
## <All variables> 0.4241204
## + Magnesium 0.4112957
## + Zinc 0.4085774
## + Nitrate_Nitrite 0.4069902
## + pH 0.4051098
## + Nickel 0.4019619
## + Phosphorus 0.4015875
## + Copper 0.3991158
## + Moisture_Content 0.3920267
## + Ammonia 0.3894283
## <none> 0.3885315
## + Calcium 0.3869749
## + Cobalt 0.3831228
##
## Df AIC F Pr(>F)
## + Magnesium 1 142.96 2.6569 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.4112957
## Call: abund_distnifH ~ Organic_Matter + Barium + Av_Potassium + Manganese + Av_Phosphorus + Iron + Sulfur + Magnesium
##
## R2.adjusted
## + Phosphorus 0.4318158
## + Nitrate_Nitrite 0.4268994
## <All variables> 0.4241204
## + Nickel 0.4189326
## + pH 0.4174772
## + Zinc 0.4166306
## + Ammonia 0.4127613
## <none> 0.4112957
## + Moisture_Content 0.4108080
## + Calcium 0.4093634
## + Copper 0.4068650
## + Cobalt 0.4064199Constrained Ordination
“http://deneflab.github.io/MicrobeMiseq/demos/mothur_2_phyloseq.html#constrained_ordinations”
CCA with selected model from Chemical CCA model
# CCA ordinate
cca_ord_nifH <- ordinate(
physeq = physeq_nifH_ord,
method = "CCA",
distance = abund_distnifH,
formula = ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium +
Av_Phosphorus + Manganese + Barium + Sulfur + Magnesium )
# CCA plot
cca_plot_nifH <- plot_ordination(
physeq = physeq_nifH_ord,
ordination = cca_ord_nifH,
color = "Site",
axes = c(1,2)) +
geom_point(aes(colour = Site), size = 3) +
scale_color_manual(values = farm_col_paired)
# Now add the environmental variables as arrows
arrowmat_nifH_cca <- vegan::scores(cca_ord_nifH, display = "bp")
#Get appropiate scalling multipler
mul<-vegan::ordiArrowMul(arrowmat_nifH_cca)
#Multiply biplot by scaling multiplier
arrowmat_nifH_cca_scale<-arrowmat_nifH_cca*3
# Add labels, make a data.frame
arrowdf_nifH_cca <- data.frame(labels = rownames(arrowmat_nifH_cca_scale), arrowmat_nifH_cca_scale)
# Define the arrow aesthetic mapping
arrow_map <- aes(xend = CCA1,
yend = CCA2,
x = 0,
y = 0,
shape = NULL,
color = NULL,
label = labels)
label_map <- aes(x = 1.3* CCA1,
y = 1.3 * CCA2,
shape = NULL,
color = NULL,
label = labels)
arrowhead = arrow(length = unit(0.02, "npc"))
# Make a new graphic
cca_plot_nifH +
geom_segment(
mapping = arrow_map,
size = .5,
data = arrowdf_nifH_cca,
color = "gray",
arrow = arrowhead
) +
geom_text(
mapping = label_map,
size = 4,
data = arrowdf_nifH_cca,
show.legend = FALSE
)+
ggtitle("CCA plot constrained ordination of nifH with selected Chemistry Model")+
theme_bw()## Warning: Ignoring unknown aesthetics: label
RDA with selected model from Chemical RDA model
# RDA ordinate
rda_ord_nifH <- ordinate(
physeq = physeq_nifH_ord,
method = "RDA",
distance = abund_distnifH,
formula = ~ Av_Phosphorus + Sulfur + Magnesium + Zinc + Av_Potassium + Manganese )
# RDA plot
rda_plot_nifH <- plot_ordination(
physeq = physeq_nifH_ord,
ordination = rda_ord_nifH,
color = "Site",
axes = c(1,2)) +
geom_point(aes(colour = Site), size = 3) +
scale_color_manual(values = farm_col_paired)
# Now add the environmental variables as arrows
arrowmat_nifH_rda <- vegan::scores(rda_ord_nifH, display = "bp")
#Get appropiate scalling multipler
mul<-vegan::ordiArrowMul(arrowmat_nifH_rda)
#Multiply biplot by scaling multiplier
arrowmat_nifH_rda_scale<-arrowmat_nifH_rda*700
# Add labels, make a data.frame
arrowdf_nifH_rda <- data.frame(labels = rownames(arrowmat_nifH_rda_scale), arrowmat_nifH_rda_scale)
# Define the arrow aesthetic mapping
arrow_map <- aes(xend = RDA1,
yend = RDA2,
x = 0,
y = 0,
shape = NULL,
color = NULL,
label = labels)
label_map <- aes(x = 1.3* RDA1,
y = 1.3 * RDA2,
shape = NULL,
color = NULL,
label = labels)
arrowhead = arrow(length = unit(0.02, "npc"))
# Make a new graphic
rda_plot_nifH +
geom_segment(
mapping = arrow_map,
size = .5,
data = arrowdf_nifH_rda,
color = "gray",
arrow = arrowhead
) +
geom_text(
mapping = label_map,
size = 4,
data = arrowdf_nifH_rda,
show.legend = FALSE
)+
ggtitle("RDA plot constrained ordination of nifH with selected Chemistry Model")+
theme_bw()## Warning: Ignoring unknown aesthetics: label
# CAP ordinate
cap_ord_nifH <- ordinate(
physeq = physeq_nifH_ord,
method = "CAP",
distance = abund_distnifH,
formula = ~ Organic_Matter + Barium + Av_Potassium +
Manganese + Av_Phosphorus + Iron + Sulfur + Magnesium )
# CAP plot
cap_plot_nifH <- plot_ordination(
physeq = physeq_nifH_ord,
ordination = cap_ord_nifH,
color = "Site",
axes = c(1,2)) +
geom_point(aes(colour = Site), size = 3) +
scale_color_manual(values = farm_col_paired)
# Now add the environmental variables as arrows
arrowmat_nifH_cap <- vegan::scores(cap_ord_nifH, display = "bp")
#Get appropiate scalling multipler
mul<-vegan::ordiArrowMul(arrowmat_nifH_cap)
#Multiply biplot by scaling multiplier
arrowmat_nifH_cap_scale<-arrowmat_nifH_cap*1.9
# Add labels, make a data.frame
arrowdf_nifH_cap <- data.frame(labels = rownames(arrowmat_nifH_cap_scale), arrowmat_nifH_cap_scale)
# Define the arrow aesthetic mapping
arrow_map <- aes(xend = CAP1,
yend = CAP2,
x = 0,
y = 0,
shape = NULL,
color = NULL,
label = labels)
label_map <- aes(x = 1.3 * CAP1,
y = 1.3 * CAP2,
shape = NULL,
color = NULL,
label = labels)
arrowhead = arrow(length = unit(0.02, "npc"))
# Make a new graphic
cap_plot_nifH +
geom_segment(
mapping = arrow_map,
size = .5,
data = arrowdf_nifH_cap,
color = "gray",
arrow = arrowhead
) +
geom_text(
mapping = label_map,
size = 4,
data = arrowdf_nifH_cap,
show.legend = FALSE
)+
ggtitle("CAP Plot Constrained Ordination of nifH with Selected Model")+
theme_bw()## Warning: Ignoring unknown aesthetics: label